#### R code for
#### Claassen & Magalhaes, 
#### Public Support for Democracy in the United States Has Declined Generationally
#### POQ
#### January 2023

#### Note: this code takes several hours to run and requires a CPU with at least 12 cores (physical or virtual)

library(arm) 
library(brms) 
library(viridis) 
library(RColorBrewer) 
library(bayesplot) 

# Options
options(mc.cores = parallel::detectCores())
par(mfrow=c(1,1), mar=c(3,3,1,1), tcl=-0.2, las=0, mgp=c(1.5, 0.4, 0))
bayesplot_theme_set(theme_default(base_size=10, base_family="sans"))
color_scheme_set(scheme="viridis")

# Folders
WD = rstudioapi::getActiveDocumentContext()$path 
setwd(dirname(WD))
print( getwd() )

# Data
# The first two datasets can be extracted from publicly available WVS and AmericasBarometer archives 
# using the "POQ US APC extract survey data.R" code
wvs = read.csv("wvs_us_democ.csv")
lapop = read.csv("lapop_us_democ.csv")
supdem = read.csv("mood_est_v5.csv")

summary(wvs)
summary(lapop)
summary(supdem)


### Exploratory analysis

## Reliability and validity of scales

wvs.it = c("eval_dem", "dem_best", "rej_milrule", "rej_strlead", "rej_expert")
psych::alpha(wvs[, wvs.it])
eigen(cor(wvs[, wvs.it], use="pair"))


## Time series plot of support in US

plot.col = viridis(n=1, begin=0.4, end=0.4)
plot.alp = viridis(n=1, begin=0.4, end=0.4, alpha=0.3)
plot.oth = RColorBrewer::brewer.pal(n=8, name="Set1")[-c(2,6)]
plot.oth = lapply(plot.oth, col2rgb)
plot.oth = lapply(plot.oth, function(x) x/255)
plot.oth = lapply(plot.oth, function(x) as.data.frame(t(x)))

cnt.plot = levels(factor(supdem$Country))
yr.surv = c(1995, 1999, 2002, 2004, 2006, 2008, 2009, 2010, 2011, 2012, 2014, 2015, 2017, 2019)

pdf("fig1a.pdf", height=5, width=6)
par(mfrow=c(1,1), mar=c(2.8, 3.2, 0.2, 0.2), cex=1, tcl=-0.3)
plot(x=supdem[supdem$Country=="United States of America", "Year"], 
     y=supdem[supdem$Country=="United States of America", "SupDem"], 
     type="n", ylim=c(-1, 2.2), axes=FALSE, xlab="", ylab="")
grid(lwd=0.66, col=rgb(0.5, 0.5, 0.5, 0.4))
rug(x=yr.surv, lwd=2, col=rgb(0, 0, 0, 1), ticksize=0.02)
lines(x=supdem[supdem$Country=="Canada", "Year"], y=supdem[supdem$Country=="Canada", "SupDem"],
      lwd=1.7, col=rgb(red=plot.oth[[1]], alpha=1), lty=2)
text(x=2000, y=supdem[supdem$Country=="Canada" & supdem$Year==2000, "SupDem"] + 0.1, "Canada", 
     col=rgb(red=plot.oth[[1]], alpha=1), cex=0.8)
lines(x=supdem[supdem$Country=="Germany", "Year"], y=supdem[supdem$Country=="Germany", "SupDem"],
      lwd=1.7, col=rgb(red=plot.oth[[2]], alpha=1), lty=2)
text(x=2020, y=supdem[supdem$Country=="Germany" & supdem$Year==2020, "SupDem"] + 0.1, "Germany", 
     col=rgb(red=plot.oth[[2]], alpha=1), cex=0.8, pos=2)
lines(x=supdem[supdem$Country=="France", "Year"], y=supdem[supdem$Country=="France", "SupDem"],
      lwd=1.7, col=rgb(red=plot.oth[[3]], alpha=1), lty=2)
text(x=1994, y=supdem[supdem$Country=="France" & supdem$Year==1995, "SupDem"] + 0.1, "France", 
     col=rgb(red=plot.oth[[3]], alpha=1), cex=0.8, pos=4)
lines(x=supdem[supdem$Country=="United Kingdom", "Year"], y=supdem[supdem$Country=="United Kingdom", "SupDem"],
      lwd=1.7, col=rgb(red=plot.oth[[4]], alpha=1), lty=2)
text(x=1994, y=supdem[supdem$Country=="United Kingdom" & supdem$Year==1995, "SupDem"] - 0.1, "UK", 
     col=rgb(red=plot.oth[[4]], alpha=1), cex=0.8, pos=4)
lines(x=supdem[supdem$Country=="Italy", "Year"], y=supdem[supdem$Country=="Italy", "SupDem"],
      lwd=1.7, col=rgb(red=plot.oth[[5]], alpha=1), lty=2)
text(x=1994, y=supdem[supdem$Country=="Italy" & supdem$Year==1995, "SupDem"] - 0.1, "Italy", 
     col=rgb(red=plot.oth[[5]], alpha=1), cex=0.8, pos=4)
lines(x=supdem[supdem$Country=="Japan", "Year"], y=supdem[supdem$Country=="Japan", "SupDem"],
      lwd=1.7, col=rgb(red=plot.oth[[6]], alpha=1), lty=2)
text(x=1996, y=supdem[supdem$Country=="Japan" & supdem$Year==1995, "SupDem"] - 0.1, "Japan", 
     col=rgb(red=plot.oth[[6]], alpha=1), cex=0.8)
lines(x=supdem[supdem$Country=="United States of America", "Year"], 
      y=supdem[supdem$Country=="United States of America", "SupDem"], 
      lwd=2.5, col=plot.col)
polygon(x=c(1995:2020, 2020:1995), 
        y=c(supdem[supdem$Country=="United States of America", "SupDem_u95"], 
            supdem[supdem$Country=="United States of America", "SupDem_l95"][26:1]),
        col=plot.alp, border=NA)
text(x=1994.5, y=supdem[supdem$Country=="United States of America" & supdem$Year==1995, "SupDem"] + 0.1, "United States", 
     col=plot.col, cex=0.9, pos=4)
axis(side=1, at=seq(1995, 2020, by=5), labels=seq(1995, 2020, by=5), 
     mgp=c(1, 0.3, 0), cex.axis=0.95)
axis(side=2, at=c(-2, -1, 0, 1, 2), mgp=c(1, 0.5, 0), cex.axis=0.95, las=1)
mtext(side=1, line=1.35, text="Year", cex=1.2)
mtext(side=2, line=1.45, text="Public support for democracy", cex=1.2)
box()
dev.off()


# WVS plot - strong leader

by(wvs$sup_dem, wvs$year, mean)

wvs$rej_strlead_di = ifelse(wvs$rej_strlead > 2, 1, 0)
summary(wvs.yr.eff <- lm(rej_strlead_di ~ -1 + factor(year), weights=wt, data=wvs))
wvs.pe = coef(wvs.yr.eff)
wvs.se = sqrt(diag(vcov(wvs.yr.eff)))
wvs.u95 = wvs.pe + 1.96*wvs.se
wvs.l95 = wvs.pe - 1.96*wvs.se
wvs.year = c(1995, 1999, 2006, 2011, 2017)

png("fig1b.png", height=400, width=600, units="px")
par(mfrow=c(1,1), mar=c(3, 3.4, 0.2, 0.2), cex=1.5, tcl=-0.4)
plot(x=wvs.year, y=wvs.pe, type="n", ylim=c(0.5, 1), xlim=c(1995, 2020), axes=FALSE, 
     xlab="", ylab="")
grid(lwd=0.66, col=rgb(0.5, 0.5, 0.5, 0.4))
points(x=wvs.year, y=wvs.pe, pch=19, cex=1.57, lwd=2, type="b", col=plot.col)
segments(x0=wvs.year, y0=wvs.u95, y1=wvs.l95, lwd=3, col=plot.col)
text(x=2020, y=0.96, labels="WVS - reject\nstrong leader", cex=1.6, adj=1)
axis(side=1, at=seq(1995, 2020, by=5), labels=seq(1995, 2020, by=5), 
     mgp=c(1, 0.4, 0), cex.axis=1.3)
axis(side=2, at=c(0.5, 0.6, 0.7, 0.8, 0.9, 1), mgp=c(1, 0.6, 0), cex.axis=1.3, las=1)
mtext(side=1, line=1.7, text="Year", cex=2.2)
mtext(side=2, line=2.3, text="Prop. rejecting strong leader", cex=2.2)
box()
dev.off()


# WVS plot appendix - strong leader

wvs$rej_strlead_di = ifelse(wvs$rej_strlead > 2, 1, 0)
summary(wvs.yr.eff <- lm(rej_strlead_di ~ -1 + factor(year), weights=wt, data=wvs))
wvs.pe = coef(wvs.yr.eff)
wvs.se = sqrt(diag(vcov(wvs.yr.eff)))
wvs.u95 = wvs.pe + 1.96*wvs.se
wvs.l95 = wvs.pe - 1.96*wvs.se
wvs.year = c(1995, 1999, 2006, 2011, 2017)

png("fig_s1a.png", height=600, width=600, units="px")
par(mfrow=c(1,1), mar=c(2.8, 3.5, 0.2, 0.2), cex=1.5, tcl=-0.4)
plot(x=wvs.year, y=wvs.pe, type="n", ylim=c(0, 1), xlim=c(1995, 2020), axes=FALSE, 
     xlab="", ylab="", yaxs="i")
grid(lwd=0.66, col=rgb(0.5, 0.5, 0.5, 0.4))
points(x=wvs.year, y=wvs.pe, pch=19, cex=1.5, type="b", col=plot.col)
segments(x0=wvs.year, y0=wvs.u95, y1=wvs.l95, lwd=2, col=plot.col)
text(x=2020, y=0.92, labels="WVS - reject\nstrong leader", cex=1.5, adj=1)
axis(side=1, at=seq(1995, 2020, by=5), labels=seq(1995, 2020, by=5), 
     mgp=c(1, 0.4, 0), cex.axis=1.2)
axis(side=2, at=c(0, 0.2, 0.4, 0.6, 0.8, 1), mgp=c(1, 0.6, 0), cex.axis=1.2, las=1)
mtext(side=1, line=1.7, text="Year", cex=2)
mtext(side=2, line=2.5, text="Proprtion rejecting strong leader", cex=2)
box()
dev.off()


# WVS plot - military rule

wvs$rej_milrule_di = ifelse(wvs$rej_milrule > 2, 1, 0)
summary(wvs.yr.eff <- lm(rej_milrule_di ~ -1 + factor(year), weights=wt, data=wvs))
wvs.pe = coef(wvs.yr.eff)
wvs.se = sqrt(diag(vcov(wvs.yr.eff)))
wvs.u95 = wvs.pe + 1.96*wvs.se
wvs.l95 = wvs.pe - 1.96*wvs.se
wvs.year = c(1995, 1999, 2006, 2011, 2017)

png("fig_s1b.png", height=600, width=600, units="px")
par(mfrow=c(1,1), mar=c(2.8, 3.5, 0.2, 0.2), cex=1.5, tcl=-0.4)
plot(x=wvs.year, y=wvs.pe, type="n", ylim=c(0, 1), xlim=c(1995, 2020), axes=FALSE, 
     xlab="", ylab="", yaxs="i")
grid(lwd=0.66, col=rgb(0.5, 0.5, 0.5, 0.4))
points(x=wvs.year, y=wvs.pe, pch=19, cex=1.5, type="b", lwd=2, col=plot.col)
segments(x0=wvs.year, y0=wvs.u95, y1=wvs.l95, col=plot.col)
text(x=2020, y=0.92, labels="WVS - reject\nmilitary rule", cex=1.5, adj=1)
axis(side=1, at=seq(1995, 2020, by=5), labels=seq(1995, 2020, by=5), 
     mgp=c(1, 0.4, 0), cex.axis=1.2)
axis(side=2, at=c(0, 0.2, 0.4, 0.6, 0.8, 1), mgp=c(1, 0.6, 0), cex.axis=1.2, las=1)
mtext(side=1, line=1.7, text="Year", cex=2)
mtext(side=2, line=2.5, text="Proprtion rejecting military rule", cex=2)
box()
dev.off()


# WVS plot - expert rule

wvs$rej_expert_di = ifelse(wvs$rej_expert > 2, 1, 0)
summary(wvs.yr.eff <- lm(rej_expert_di ~ -1 + factor(year), weights=wt, data=wvs))
wvs.pe = coef(wvs.yr.eff)
wvs.se = sqrt(diag(vcov(wvs.yr.eff)))
wvs.u95 = wvs.pe + 1.96*wvs.se
wvs.l95 = wvs.pe - 1.96*wvs.se
wvs.year = c(1995, 1999, 2006, 2011, 2017)

png("fig_s1c.png", height=600, width=600, units="px")
par(mfrow=c(1,1), mar=c(2.8, 3.5, 0.2, 0.2), cex=1.5, tcl=-0.4)
plot(x=wvs.year, y=wvs.pe, type="n", ylim=c(0, 1), xlim=c(1995, 2020), axes=FALSE, 
     xlab="", ylab="", yaxs="i")
grid(lwd=0.66, col=rgb(0.5, 0.5, 0.5, 0.4))
points(x=wvs.year, y=wvs.pe, pch=19, cex=1.5, type="b", col=plot.col)
segments(x0=wvs.year, y0=wvs.u95, y1=wvs.l95, lwd=2, col=plot.col)
text(x=2020, y=0.92, labels="WVS - reject\nexpert rule", cex=1.5, adj=1)
axis(side=1, at=seq(1995, 2020, by=5), labels=seq(1995, 2020, by=5), 
     mgp=c(1, 0.4, 0), cex.axis=1.2)
axis(side=2, at=c(0, 0.2, 0.4, 0.6, 0.8, 1), mgp=c(1, 0.6, 0), cex.axis=1.2, las=1)
mtext(side=1, line=1.7, text="Year", cex=2)
mtext(side=2, line=2.5, text="Proprtion rejecting expert rule", cex=2)
box()
dev.off()


# WVS plot - support democracy

wvs$eval_dem_di = ifelse(wvs$eval_dem > 2, 1, 0)
summary(wvs.yr.eff <- lm(eval_dem_di ~ -1 + factor(year), weights=wt, data=wvs))
wvs.pe = coef(wvs.yr.eff)
wvs.se = sqrt(diag(vcov(wvs.yr.eff)))
wvs.u95 = wvs.pe + 1.96*wvs.se
wvs.l95 = wvs.pe - 1.96*wvs.se
wvs.year = c(1995, 1999, 2006, 2011, 2017)

png("fig_s1d.png", height=600, width=600, units="px")
par(mfrow=c(1,1), mar=c(2.8, 3.5, 0.2, 0.2), cex=1.5, tcl=-0.4)
plot(x=wvs.year, y=wvs.pe, type="n", ylim=c(0, 1), xlim=c(1995, 2020), axes=FALSE, 
     xlab="", ylab="", yaxs="i")
grid(lwd=0.66, col=rgb(0.5, 0.5, 0.5, 0.4))
points(x=wvs.year, y=wvs.pe, pch=19, cex=1.5, type="b", col=plot.col)
segments(x0=wvs.year, y0=wvs.u95, y1=wvs.l95, lwd=2, col=plot.col)
text(x=2020, y=0.92, labels="WVS - support\ndemocracy", cex=1.5, adj=1)
axis(side=1, at=seq(1995, 2020, by=5), labels=seq(1995, 2020, by=5), 
     mgp=c(1, 0.4, 0), cex.axis=1.2)
axis(side=2, at=c(0, 0.2, 0.4, 0.6, 0.8, 1), mgp=c(1, 0.6, 0), cex.axis=1.2, las=1)
mtext(side=1, line=1.7, text="Year", cex=2)
mtext(side=2, line=2.5, text="Proprtion rating democracy good", cex=2)
box()
dev.off()


# AB plot

by(lapop$dem_best, lapop$year, table)

summary(ab.yr.eff <- lm(dem_best ~ -1 + factor(year), weights=wt, data=lapop))
ab.pe = coef(ab.yr.eff)
ab.se = sqrt(diag(vcov(ab.yr.eff)))
ab.u95 = ab.pe + 1.96*ab.se
ab.l95 = ab.pe - 1.96*ab.se
ab.year = c(2006, 2008, 2010, 2012, 2014, 2017, 2019)

lapop$dem_best_di = ifelse(lapop$dem_best > 4, 1, 0)
summary(ab.yr.eff <- lm(dem_best_di ~ -1 + factor(year), weights=wt, data=lapop))
ab.pe = coef(ab.yr.eff)
ab.se = sqrt(diag(vcov(ab.yr.eff)))
ab.u95 = ab.pe + 1.96*ab.se
ab.l95 = ab.pe - 1.96*ab.se
ab.year = c(2006, 2008, 2010, 2012, 2014, 2017, 2019)

png("fig1c.png", height=400, width=600, units="px")
par(mfrow=c(1,1), mar=c(3, 3.4, 0.2, 0.2), cex=1.5, tcl=-0.4)
plot(x=ab.year, y=ab.pe, type="n", ylim=c(0.5, 1), xlim=c(1995, 2020), axes=FALSE, 
     xlab="", ylab="")
grid(lwd=0.66, col=rgb(0.5, 0.5, 0.5, 0.4))
points(x=ab.year, y=ab.pe, pch=19, cex=1.57, type="b", lwd=2, col=plot.col)
segments(x0=ab.year, y0=ab.u95, y1=ab.l95, lwd=3, col=plot.col)
text(x=2020, y=0.96, labels="AmericasBarometer\n- democracy is best", cex=1.6, adj=1)
axis(side=1, at=seq(1995, 2020, by=5), labels=seq(1995, 2020, by=5), 
     mgp=c(1, 0.4, 0), cex.axis=1.3)
axis(side=2, at=c(0.5, 0.6, 0.7, 0.8, 0.9, 1), mgp=c(1, 0.6, 0), cex.axis=1.3, las=1)
mtext(side=1, line=1.7, text="Year", cex=2.2)
mtext(side=2, line=2.3, text="Prop. supporting democracy", cex=2.2)
box()
dev.off()


## Bivariate plots of generations

# WVS

dec_mod = lm(sup_dem ~ -1 + as.factor(decbirth), wvs, weights=wt)
wvs_gen_se = sqrt(diag(vcov(dec_mod)))
wvs_gen_pe = coef(dec_mod)
wvs_gen_up = wvs_gen_pe + 1.96*wvs_gen_se
wvs_gen_lo = wvs_gen_pe - 1.96*wvs_gen_se

plot.col = inferno(n=1, begin=0.7, end=0.7)

pdf("fig_s2b.pdf", height=4, width=5.5)
wvs_gen = c("1900s","1910s","1920s","1930s","1940s","1950s","1960s","1970s","1980s","1990s")
par(mfrow=c(1,1), mar=c(3, 3, 2.2, 0.5), tcl=-0.2, las=0, cex.axis=0.9)
plot(1:10, y=wvs_gen_pe, type="n", ylim=c(-1, 1.1), xlab="", ylab="", axes=FALSE)
abline(h=seq(-2, 1.5, by=0.5), col=rgb(0.5, 0.5, 0.5, 0.9), lwd=0.75, lty=3)
abline(v=1:10, col=rgb(0.5, 0.5, 0.5, 0.9), lwd=0.75, lty=3)
segments(x0=1:10, y0=wvs_gen_up, y1=wvs_gen_lo, lwd=2, col=plot.col)
points(x=1:10, y=wvs_gen_pe, col=plot.col, pch=16, cex=1.2)
axis(side=1, at=1:10, labels=wvs_gen, cex.axis=0.99, mgp=c(1, 0.15, 0))
axis(side=2, at=seq(-1, 1, by=0.5), cex.axis=0.99, las=1, mgp=c(1, 0.35, 0))
text(x=9.8, y=1.05, "WVS", cex=0.9)
mtext(side=1, line=1.3, "Decade of birth", cex=0.9)
mtext(side=2, line=1.5, text="Support for democracy (5-item), std scale", cex=1, las=0)
mtext(side=3, line=0.5, text="Support for democracy by birth decade, bivariate", cex=1.2, las=0)
box()
dev.off()


# LAPOP

l.g.prior = c(set_prior("normal(0, 1)", class="b"),
              set_prior("normal(0, 2)", class="Intercept"))

lapop2 = lapop
lapop2$dem_best = as.ordered(lapop2$dem_best)
lapop2$birthdec = as.factor(lapop2$birthdec)
l.g.ser <- brm(dem_best_ord | weights(wt) ~ birthdec, data=lapop2, family=cumulative("logit"),
               prior=l.g.prior, control=list(adapt_delta=0.99), chains=12, iter=12, backend="cmdstanr")

l.g.par <- update(
  l.g.ser, chains=3, cores=3, iter=2000, backend="cmdstanr", threads=threading(4),
  control=c(adapt_delta=0.99, max_treedepth=12), stan_model_args=list(stanc_options = list("O1"))
)

summary(l.g.par)
pp_check(l.g.par)

gen_labs = c("1910s", "1920s", "1930s", "1940s", "1950s", "1960s", "1970s", "1980s", "1990s", "2000s")

pred.dat.l.g = data.frame(birthdec=1:10)
pred.prob.l.g = predict(l.g.par, newdata=pred.dat.l.g, cores=8, summary=TRUE)
rowSums(pred.prob.l.g[,5:7])

pred.arr.l.g = posterior_epred(l.g.par, newdata=pred.dat.l.g, cores=8)
pred.arr.l.g = pred.arr.l.g[,,5:7]
pred.l.g = apply(pred.arr.l.g, c(1,2), sum)
gen_pe = apply(pred.l.g, 2, mean)
gen_lo = apply(pred.l.g, 2, function(x) quantile(x, probs=0.025))
gen_up = apply(pred.l.g, 2, function(x) quantile(x, probs=0.975))

plot.col = inferno(n=1, begin=0.7, end=0.7)
plot.alp = inferno(n=1, begin=0.7, end=0.7, alpha=0.5)
 
pdf("fig_s2a.pdf", height=4, width=5.5)
par(mfrow=c(1,1), mar=c(3, 3, 2.2, 0.5), tcl=-0.2, las=0)
plot(x=1:10, y=gen_pe, pch=16, xlim=c(1, 10), ylim=c(0.4, 1), type="n",
     xlab="", ylab="", main="", axes=FALSE, yaxs="i")
abline(h=seq(0.4, 1, by=0.1), col=rgb(0.5, 0.5, 0.5, 0.9), lwd=0.75, lty=3)
abline(v=1:10, col=rgb(0.5, 0.5, 0.5, 0.9), lwd=0.75, lty=3)
segments(x0=1:10, y0=gen_up, y1=gen_lo, lwd=2, col=plot.col)
points(x=1:10, y=gen_pe, col=plot.col, pch=16, cex=1.2)
axis(1, at=1:10, labels=gen_labs, cex.axis=0.99, las=1, mgp=c(1, 0.15, 0))
axis(2, at=seq(0.4, 1, by=0.1), cex.axis=0.99, las=1, mgp=c(1, 0.35, 0))
text(x=10, y=0.98, labels="AmericasBarometer", pos=2, cex=1)
mtext(side=1, line=1.3, "Decade of birth", cex=1.1)
mtext(side=2, line=1.5, "Probability of agreeing democracy is best", cex=1)
mtext("Support for democracy by birth decade, bivariate", side=3, line=0.5, adj=0.5, cex=1.2)
box()
dev.off()



#### Age Period Cohort models


### LAPOP

## Hierarchical ordinal APC models - Bayesian MCMC - 4 age groups - w/ controls - in SI

h.b.l.c.prior = c(set_prior("normal(0, 1)", class="b"),
                  set_prior("normal(0, 2)", class="Intercept"),
                  set_prior("normal(0, 1)", class="sd"))

lapop$dem_best_ord = as.ordered(lapop$dem_best)
h.b.l.c.1.ser <- brm(dem_best_ord | weights(wt) ~ republican + democrat + degree + female + white + income_std 
                     + (1 | year) + (1 | agecat4) + (1 | birthdec) + south, data=lapop, family=cumulative("logit"), 
                     prior=h.b.l.c.prior, control=list(adapt_delta=0.99), chains=12, iter=12, backend="cmdstanr")

h.b.l.c.1.par <- update(
  h.b.l.c.1.ser, chains=3, cores=3, iter=2000, backend="cmdstanr", threads=threading(4),
  control=c(adapt_delta=0.99, max_treedepth=12), stan_model_args=list(stanc_options = list("O1"))
)

summary(h.b.l.c.1.par)
pp_check(h.b.l.c.1.par)

pdf("fig_s9b.pdf", width=6, height=6)
plot(h.b.l.c.1.par, newpage=FALSE, 
     variable=c("b_Intercept[1]", "sd_agecat4__Intercept", "sd_birthdec__Intercept", "sd_year__Intercept"))
dev.off()

sum.h.b.l.c.1.par <- capture.output(summary(h.b.l.c.1.par))
cat(title="Lapop HAPCM results", sum.h.b.l.c.1.par, file="table_s2_2.txt", sep="\n", append=TRUE)

plot(ranef(h.b.l.c.1.par)$year[,,1][,1], type="b", xlab="Year", ylab="Support democracy")
plot(ranef(h.b.l.c.1.par)$agecat4[,,1][,1], type="b", xlab="Age category", ylab="Support democracy")
plot(ranef(h.b.l.c.1.par)$birthdec[,,1][,1], type="b", xlab="Birth decade", ylab="Support democracy")

# cohort plots

gen_pe = 1- invlogit(coef(h.b.l.c.1.par)$birthdec[,,4][,1])
gen_up = 1- invlogit(coef(h.b.l.c.1.par)$birthdec[,,4][,4])
gen_lo = 1- invlogit(coef(h.b.l.c.1.par)$birthdec[,,4][,3])
gen_labs = c("1910s", "1920s", "1930s", "1940s", "1950s", "1960s", "1970s", "1980s", "1990s", "2000s")

pred.dat.h.b.l.c.1 = data.frame(republican=0, democrat=0, degree=0, female=1, white=1, income_std=0, year=2019,
                                agecat4=2, birthdec=1:10, south=0)
pred.prob.h.b.l.c.1 = predict(h.b.l.c.1.par, newdata=pred.dat.h.b.l.c.1, cores=8, summary=TRUE)
rowSums(pred.prob.h.b.l.c.1[,5:7])

pred.arr.h.b.l.c.1 = posterior_epred(h.b.l.c.1.par, newdata=pred.dat.h.b.l.c.1, cores=8)
pred.arr.h.b.l.c.1 = pred.arr.h.b.l.c.1[,,5:7]
pred.arr.h.b.l.c.1 = apply(pred.arr.h.b.l.c.1, c(1,2), sum)
gen_pe = apply(pred.arr.h.b.l.c.1, 2, mean)
gen_lo = apply(pred.arr.h.b.l.c.1, 2, function(x) quantile(x, probs=0.025))
gen_up = apply(pred.arr.h.b.l.c.1, 2, function(x) quantile(x, probs=0.975))

pdf("fig_s3c.pdf", height=4, width=5.5)
plot.col = viridis(n=1, begin=0.4, end=0.4)
plot.alp = viridis(n=1, begin=0.4, end=0.4, alpha=0.5)
par(mfrow=c(1,1), mar=c(3, 3, 2.2, 0.5), tcl=-0.2, las=0)
plot(x=1:10, y=gen_pe, pch=16, xlim=c(1, 10), ylim=c(0.3, 1), type="n",
     xlab="", ylab="", main="", axes=FALSE, yaxs="i")
abline(h=seq(0.4, 1, by=0.1), col=rgb(0.5, 0.5, 0.5, 0.9), lwd=0.75, lty=3)
abline(v=1:10, col=rgb(0.5, 0.5, 0.5, 0.9), lwd=0.75, lty=3)
lines(x=1:10, y=gen_pe, lwd=2, col=plot.col)
polygon(x=c(1:10, 10:1), y=c(gen_up, gen_lo[10:1]), col=plot.alp, border=NA, density=NA)
axis(1, at=1:10, labels=gen_labs, cex.axis=0.99, las=1, mgp=c(1, 0.15, 0))
axis(2, at=seq(0.4, 1, by=0.1), cex.axis=0.99, las=1, mgp=c(1, 0.35, 0))
text(x=10, y=0.98, labels="AmericasBarometer", pos=2, cex=1)
mtext(side=1, line=1.3, "Decade of birth", cex=1.1)
mtext(side=2, line=1.5, "Probability of agreeing democracy is best", cex=1)
mtext("Support for democracy by birth decade, controls", side=3, line=0.5, adj=0.5, cex=1.2)
box()
dev.off()

# age plots

age_pe = ranef(h.b.l.c.1.par)$agecat4[,,1][,1]
age_up = ranef(h.b.l.c.1.par)$agecat4[,,1][,4]
age_lo = ranef(h.b.l.c.1.par)$agecat4[,,1][,3]
age_labs = c("18-29", "30-44", "45-59", "60+")

pred.dat.h.b.l.c.1 = data.frame(year=2019, agecat4=1:4, birthdec=6, republican=0, democrat=0, degree=0, 
                                female=1, white=1, income_std=0, south=0)
pred.prob.h.b.l.c.1 = predict(h.b.l.c.1.par, newdata=pred.dat.h.b.l.c.1, cores=8, summary=TRUE)
rowSums(pred.prob.h.b.l.c.1[,5:7])

pred.arr.h.b.l.c.1 = posterior_epred(h.b.l.c.1.par, newdata=pred.dat.h.b.l.c.1, cores=8)
pred.arr.h.b.l.c.1 = pred.arr.h.b.l.c.1[,,5:7]
pred.arr.h.b.l.c.1 = apply(pred.arr.h.b.l.c.1, c(1,2), sum)
age_pe = apply(pred.arr.h.b.l.c.1, 2, mean)
age_lo = apply(pred.arr.h.b.l.c.1, 2, function(x) quantile(x, probs=0.025))
age_up = apply(pred.arr.h.b.l.c.1, 2, function(x) quantile(x, probs=0.975))

pdf("fig_s4c.pdf", height=4, width=5.5)
plot.col = viridis(n=1, begin=0.4, end=0.4)
plot.alp = viridis(n=1, begin=0.4, end=0.4, alpha=0.5)
par(mfrow=c(1,1), mar=c(3, 3, 2.2, 0.5), tcl=-0.2, las=0)
plot(x=1:4, y=age_pe, pch=16, xlim=c(0.8, 4.2), ylim=c(0.55, 0.85), type="n",
     xlab="", ylab="", main="", axes=FALSE)
abline(h=seq(0.5, 0.8, by=0.1), col=rgb(0.5, 0.5, 0.5, 0.9), lwd=0.75, lty=3)
abline(v=1:4, col=rgb(0.5, 0.5, 0.5, 0.9), lwd=0.75, lty=3)
points(x=1:4, y=age_pe, pch=16, cex=1.2, col=plot.col)
segments(x0=1:4, y0=age_lo, y1=age_up, lwd=2, col=plot.col)
axis(1, at=1:4, labels=age_labs, cex.axis=0.99, las=1, mgp=c(1, 0.15, 0))
axis(2, at=seq(0.5, 0.8, by=0.1), cex.axis=0.99, las=1, mgp=c(1, 0.35, 0))
text(x=4.2, y=0.85, labels="AmericasBarometer", pos=2, cex=1)
mtext(side=1, line=1.3, "Age", cex=1.1)
mtext(side=2, line=1.5, "Probability of agreeing democracy is best", cex=1)
mtext("Support for democracy by age, controls", side=3, line=0.5, adj=0.5, cex=1.2)
box()
dev.off()

# period plots

year_pe = ranef(h.b.l.c.1.par)$year[,,1][,1]
year_up = ranef(h.b.l.c.1.par)$year[,,1][,4]
year_lo = ranef(h.b.l.c.1.par)$year[,,1][,3]
year_labs = c("2006", "2008", "2010", "2012", "2014", "2017", "2019")

pred.dat.h.b.l.c.1 = data.frame(year=as.numeric(year_labs), agecat4=3, birthdec=6, republican=0, 
                                democrat=0, degree=0, female=1, white=1, income_std=0, south=0)
pred.prob.h.b.l.c.1 = predict(h.b.l.c.1.par, newdata=pred.dat.h.b.l.c.1, cores=8, summary=TRUE)
rowSums(pred.prob.h.b.l.c.1[,5:7])

pred.arr.h.b.l.c.1 = posterior_epred(h.b.l.c.1.par, newdata=pred.dat.h.b.l.c.1, cores=8)
pred.arr.h.b.l.c.1 = pred.arr.h.b.l.c.1[,,5:7]
pred.arr.h.b.l.c.1 = apply(pred.arr.h.b.l.c.1, c(1,2), sum)
year_pe = apply(pred.arr.h.b.l.c.1, 2, mean)
year_lo = apply(pred.arr.h.b.l.c.1, 2, function(x) quantile(x, probs=0.025))
year_up = apply(pred.arr.h.b.l.c.1, 2, function(x) quantile(x, probs=0.975))

pdf("fig_s5c.pdf", height=4, width=5.5)
plot.col = viridis(n=1, begin=0.4, end=0.4)
plot.alp = viridis(n=1, begin=0.4, end=0.4, alpha=0.5)
par(mfrow=c(1,1), mar=c(3, 3, 2.2, 0.5), tcl=-0.2, las=0)
plot(x=1:7, y=year_pe, pch=16, xlim=c(0.8, 7.2), ylim=c(0.5, 1), type="n",
     xlab="", ylab="", main="", axes=FALSE, yaxs="i")
abline(h=seq(0.5, 1, by=0.1), col=rgb(0.5, 0.5, 0.5, 0.9), lwd=0.75, lty=3)
abline(v=1:7, col=rgb(0.5, 0.5, 0.5, 0.9), lwd=0.75, lty=3)
points(x=1:7, y=year_pe, pch=16, cex=1.2, col=plot.col)
segments(x0=1:7, y0=year_lo, y1=year_up, lwd=2, col=plot.col)
axis(1, at=1:7, labels=year_labs, cex.axis=0.99, las=1, mgp=c(1, 0.15, 0))
axis(2, at=seq(0.5, 1, by=0.1), cex.axis=0.99, las=1, mgp=c(1, 0.35, 0))
text(x=7.2, y=0.98, labels="AmericasBarometer", pos=2, cex=1)
mtext(side=1, line=1.3, "Year", cex=1.1)
mtext(side=2, line=1.6, "Probability of agreeing democracy is best", cex=1)
mtext("Support for democracy by year, controls", side=3, line=0.5, adj=0.5, cex=1.2)
box()
dev.off()


## Hierarchical ordinal APC models - Bayesian MCMC - 4 age groups - no controls - in SI

h.b.l.nc.prior = c(set_prior("normal(0, 2)", class="Intercept"),
                  set_prior("normal(0, 1)", class="sd"))

lapop$dem_best_ord = as.ordered(lapop$dem_best)
h.b.l.nc.1.ser <- brm(dem_best_ord | weights(wt) ~ (1 | year) + (1 | agecat4) + (1 | birthdec), 
                      data=lapop, family=cumulative("logit"), prior=h.b.l.nc.prior, 
                      control=list(adapt_delta=0.99), chains=12, iter=12, backend="cmdstanr")

h.b.l.nc.1.par <- update(
  h.b.l.nc.1.ser, chains=3, cores=3, iter=2000, backend="cmdstanr", threads=threading(4),
  control=c(adapt_delta=0.99, max_treedepth=12), stan_model_args=list(stanc_options = list("O1"))
)

summary(h.b.l.nc.1.par)
pp_check(h.b.l.nc.1.par)

pdf("fig_s9a.pdf", width=6, height=6)
plot(h.b.l.nc.1.par, newpage=FALSE, 
     variable=c("b_Intercept[1]", "sd_agecat4__Intercept", "sd_birthdec__Intercept", "sd_year__Intercept"))
dev.off()

sum.h.b.l.nc.1.par <- capture.output(summary(h.b.l.nc.1.par))
cat(title="Lapop HAPCM results", sum.h.b.l.nc.1.par, file="table_s2_1.txt", sep="\n", append=TRUE)

plot(ranef(h.b.l.nc.1.par)$year[,,1][,1], type="b", xlab="Year", ylab="Support democracy")
plot(ranef(h.b.l.nc.1.par)$agecat4[,,1][,1], type="b", xlab="Age category", ylab="Support democracy")
plot(ranef(h.b.l.nc.1.par)$birthdec[,,1][,1], type="b", xlab="Birth decade", ylab="Support democracy")

# cohort plots

gen_pe = 1- invlogit(coef(h.b.l.nc.1.par)$birthdec[,,4][,1])
gen_up = 1- invlogit(coef(h.b.l.nc.1.par)$birthdec[,,4][,4])
gen_lo = 1- invlogit(coef(h.b.l.nc.1.par)$birthdec[,,4][,3])
gen_labs = c("1910s", "1920s", "1930s", "1940s", "1950s", "1960s", "1970s", "1980s", "1990s", "2000s")

pred.dat.h.b.l.nc.1 = data.frame(year=2019, agecat4=2, birthdec=1:10)
pred.prob.h.b.l.nc.1 = predict(h.b.l.nc.1.par, newdata=pred.dat.h.b.l.nc.1, cores=8, summary=TRUE)
rowSums(pred.prob.h.b.l.nc.1[,5:7])

pred.arr.h.b.l.nc.1 = posterior_epred(h.b.l.nc.1.par, newdata=pred.dat.h.b.l.nc.1, cores=8)
pred.arr.h.b.l.nc.1 = pred.arr.h.b.l.nc.1[,,5:7]
pred.arr.h.b.l.nc.1 = apply(pred.arr.h.b.l.nc.1, c(1,2), sum)
gen_pe = apply(pred.arr.h.b.l.nc.1, 2, mean)
gen_lo = apply(pred.arr.h.b.l.nc.1, 2, function(x) quantile(x, probs=0.025))
gen_up = apply(pred.arr.h.b.l.nc.1, 2, function(x) quantile(x, probs=0.975))

pdf("fig_s3a.pdf", height=4, width=5.5)
plot.col = viridis(n=1, begin=0.4, end=0.4)
plot.alp = viridis(n=1, begin=0.4, end=0.4, alpha=0.5)
par(mfrow=c(1,1), mar=c(3, 3, 2.2, 0.5), tcl=-0.2, las=0)
plot(x=1:10, y=gen_pe, pch=16, xlim=c(1, 10), ylim=c(0.3, 1), type="n",
     xlab="", ylab="", main="", axes=FALSE, yaxs="i")
abline(h=seq(0.4, 1, by=0.1), col=rgb(0.5, 0.5, 0.5, 0.9), lwd=0.75, lty=3)
abline(v=1:10, col=rgb(0.5, 0.5, 0.5, 0.9), lwd=0.75, lty=3)
lines(x=1:10, y=gen_pe, lwd=2, col=plot.col)
polygon(x=c(1:10, 10:1), y=c(gen_up, gen_lo[10:1]), col=plot.alp, border=NA, density=NA)
axis(1, at=1:10, labels=gen_labs, cex.axis=0.99, las=1, mgp=c(1, 0.15, 0))
axis(2, at=seq(0.4, 1, by=0.1), cex.axis=0.99, las=1, mgp=c(1, 0.35, 0))
text(x=10, y=0.98, labels="AmericasBarometer", pos=2, cex=1)
mtext(side=1, line=1.3, "Decade of birth", cex=1.1)
mtext(side=2, line=1.5, "Probability of agreeing democracy is best", cex=1)
mtext("Support for democracy by birth decade, no controls", side=3, line=0.5, adj=0.5, cex=1.2)
box()
dev.off()

# age plots

age_pe = ranef(h.b.l.nc.1.par)$agecat4[,,1][,1]
age_up = ranef(h.b.l.nc.1.par)$agecat4[,,1][,4]
age_lo = ranef(h.b.l.nc.1.par)$agecat4[,,1][,3]
age_labs = c("18-29", "30-44", "45-59", "60+")

pred.dat.h.b.l.nc.1 = data.frame(year=2019, agecat4=1:4, birthdec=6)
pred.prob.h.b.l.nc.1 = predict(h.b.l.nc.1.par, newdata=pred.dat.h.b.l.nc.1, cores=8, summary=TRUE)
rowSums(pred.prob.h.b.l.nc.1[,5:7])

pred.arr.h.b.l.nc.1 = posterior_epred(h.b.l.nc.1.par, newdata=pred.dat.h.b.l.nc.1, cores=8)
pred.arr.h.b.l.nc.1 = pred.arr.h.b.l.nc.1[,,5:7]
pred.arr.h.b.l.nc.1 = apply(pred.arr.h.b.l.nc.1, c(1,2), sum)
age_pe = apply(pred.arr.h.b.l.nc.1, 2, mean)
age_lo = apply(pred.arr.h.b.l.nc.1, 2, function(x) quantile(x, probs=0.025))
age_up = apply(pred.arr.h.b.l.nc.1, 2, function(x) quantile(x, probs=0.975))

pdf("fig_s4a.pdf", height=4, width=5.5)
plot.col = viridis(n=1, begin=0.4, end=0.4)
plot.alp = viridis(n=1, begin=0.4, end=0.4, alpha=0.5)
par(mfrow=c(1,1), mar=c(3, 3, 2.2, 0.5), tcl=-0.2, las=0)
plot(x=1:4, y=age_pe, pch=16, xlim=c(0.8, 4.2), ylim=c(0.55, 0.85), type="n",
     xlab="", ylab="", main="", axes=FALSE)
abline(h=seq(0.6, 0.8, by=0.1), col=rgb(0.5, 0.5, 0.5, 0.9), lwd=0.75, lty=3)
abline(v=1:4, col=rgb(0.5, 0.5, 0.5, 0.9), lwd=0.75, lty=3)
points(x=1:4, y=age_pe, pch=16, cex=1.2, col=plot.col)
segments(x0=1:4, y0=age_lo, y1=age_up, lwd=2, col=plot.col)
axis(1, at=1:4, labels=age_labs, cex.axis=0.99, las=1, mgp=c(1, 0.15, 0))
axis(2, at=seq(0.6, 0.8, by=0.1), cex.axis=0.99, las=1, mgp=c(1, 0.35, 0))
text(x=4.2, y=0.85, labels="AmericasBarometer", pos=2, cex=1)
mtext(side=1, line=1.3, "Age", cex=1.1)
mtext(side=2, line=1.5, "Probability of agreeing democracy is best", cex=1)
mtext("Support for democracy by age, no controls", side=3, line=0.5, adj=0.5, cex=1.2)
box()
dev.off()

# period plots

year_pe = ranef(h.b.l.nc.1.par)$year[,,1][,1]
year_up = ranef(h.b.l.nc.1.par)$year[,,1][,4]
year_lo = ranef(h.b.l.nc.1.par)$year[,,1][,3]
year_labs = c("2006", "2008", "2010", "2012", "2014", "2017", "2019")

pred.dat.h.b.l.nc.1 = data.frame(year=as.numeric(year_labs), agecat4=3, birthdec=6)
pred.prob.h.b.l.nc.1 = predict(h.b.l.nc.1.par, newdata=pred.dat.h.b.l.nc.1, cores=8, summary=TRUE)
rowSums(pred.prob.h.b.l.nc.1[,5:7])

pred.arr.h.b.l.nc.1 = posterior_epred(h.b.l.nc.1.par, newdata=pred.dat.h.b.l.nc.1, cores=8)
pred.arr.h.b.l.nc.1 = pred.arr.h.b.l.nc.1[,,5:7]
pred.arr.h.b.l.nc.1 = apply(pred.arr.h.b.l.nc.1, c(1,2), sum)
year_pe = apply(pred.arr.h.b.l.nc.1, 2, mean)
year_lo = apply(pred.arr.h.b.l.nc.1, 2, function(x) quantile(x, probs=0.025))
year_up = apply(pred.arr.h.b.l.nc.1, 2, function(x) quantile(x, probs=0.975))

pdf("fig_s5a.pdf", height=4, width=5.5)
plot.col = viridis(n=1, begin=0.4, end=0.4)
plot.alp = viridis(n=1, begin=0.4, end=0.4, alpha=0.5)
par(mfrow=c(1,1), mar=c(3, 3, 2.2, 0.5), tcl=-0.2, las=0)
plot(x=1:7, y=year_pe, pch=16, xlim=c(0.8, 7.2), ylim=c(0.5, 1), type="n",
     xlab="", ylab="", main="", axes=FALSE, yaxs="i")
abline(h=seq(0.5, 1, by=0.1), col=rgb(0.5, 0.5, 0.5, 0.9), lwd=0.75, lty=3)
abline(v=1:7, col=rgb(0.5, 0.5, 0.5, 0.9), lwd=0.75, lty=3)
points(x=1:7, y=year_pe, pch=16, cex=1.2, col=plot.col)
segments(x0=1:7, y0=year_lo, y1=year_up, lwd=2, col=plot.col)
axis(1, at=1:7, labels=year_labs, cex.axis=0.99, las=1, mgp=c(1, 0.15, 0))
axis(2, at=seq(0.6, 1, by=0.1), cex.axis=0.99, las=1, mgp=c(1, 0.35, 0))
text(x=7.2, y=0.98, labels="AmericasBarometer", pos=2, cex=1)
mtext(side=1, line=1.3, "Year", cex=1.1)
mtext(side=2, line=1.6, "Probability of agreeing democracy is best", cex=1)
mtext("Support for democracy by year, no controls", side=3, line=0.5, adj=0.5, cex=1.2)
box()
dev.off()


## Bayesian ordinal GAMs using BRMS - democracy is best / LAPOP - with controls - in paper

g.o.b.l.prior = c(set_prior("normal(0, 2)", class="b"),
                  set_prior("normal(0, 5)", class="Intercept"),
                  set_prior("normal(0, 1)", class="sds"))

lapop2 = lapop
lapop2$year = as.factor(lapop$year)
lapop2$agecat4 = as.factor(lapop$agecat4)
lapop2$dem_best = as.ordered(lapop$dem_best)

g.o.b.l.c.2.ser <- brm(dem_best | weights(wt) ~ s(birthyr) + republican + democrat + degree + female + white 
                       + income_std + south + year + agecat4, data=lapop2, prior=g.o.b.l.prior, 
                       family=cumulative("logit"), control=list(adapt_delta=0.95), chains=6, iter=12, 
                       backend="cmdstanr")

g.o.b.l.c.2.par <- update(
  g.o.b.l.c.2.ser, chains=3, cores=3, iter=2000, backend="cmdstanr", threads=threading(4),
  control=c(adapt_delta=0.95, max_treedepth=12), stan_model_args=list(stanc_options = list("O1"))
)

summary(g.o.b.l.c.2.par)
pp_check(g.o.b.l.c.2.par)
pp_check(g.o.b.l.c.2.par, type="ecdf_overlay")

pdf("fig_s7b.pdf", width=6, height=6)
plot(g.o.b.l.c.2.par, newpage=FALSE, 
     variable=c("b_Intercept[1]", "b_year2008", "b_agecat42", "bs_sbirthyr_1", "sds_sbirthyr_1"))
dev.off()

sum.g.o.b.l.c.2.par <- capture.output(summary(g.o.b.l.c.2.par))
cat(title="Lapop GAM results", sum.g.o.b.l.c.2.par, file="table2_2.txt", sep="\n", append=TRUE)

# age plot

pred.dat.age = data.frame(birthyr=1964, republican=0, democrat=0, degree=0, female=1, white=1, income_std=0, 
                          south=0, year=2012, agecat4=1:4)
g.o.b.l.c.2.age.pred.eprob = posterior_epred(g.o.b.l.c.2.par, newdata=pred.dat.age, ndraws=1000)
g.o.b.l.c.2.age.pred.eprob2 = g.o.b.l.c.2.age.pred.eprob[,,5] + g.o.b.l.c.2.age.pred.eprob[,,6] + g.o.b.l.c.2.age.pred.eprob[,,7]
g.o.b.l.c.2.age.pred.mean = apply(g.o.b.l.c.2.age.pred.eprob2, 2, mean)
g.o.b.l.c.2.age.pred.l95 = apply(g.o.b.l.c.2.age.pred.eprob2, 2, quantile, probs=0.025)
g.o.b.l.c.2.age.pred.u95 = apply(g.o.b.l.c.2.age.pred.eprob2, 2, quantile, probs=0.975)

pdf("fig3c.pdf", height=4, width=5.5)
par(mfrow=c(1,1), mar=c(3, 3, 2.2, 0.5), tcl=-0.2, las=0, cex=1)
plot(x=1:4, y=g.o.b.l.c.2.age.pred.mean, type="n", axes=FALSE, xlab="", ylab="", xlim=c(0.8, 4.2), 
     ylim=c(0.6, 0.85))
axis(side=1, at=1:4, labels=c("18-29", "30-44", "45-59", "60+"), mgp=c(1, 0.1, 0), cex.axis=0.99)
axis(side=2, at=seq(0.6, 0.9, by=0.1), mgp=c(1, 0.3, 0), cex.axis=0.99)
mtext(side=1, line=1.3, "Age group", cex=1.1)
mtext(side=2, line=1.4, "Probability of agreeing democracy is best (5-7)", cex=1)
mtext(side=3, line=0.5, "Support for democracy by age group, controls", cex=1.2, adj=0.5)
abline(h=seq(0.6, 0.8, by=0.05), lwd=0.75, lty=3, col=rgb(0.5, 0.5, 0.5, 0.7))
abline(h=0, lwd=0.75, col=rgb(0.5, 0.5, 0.5, 0.6))
abline(v=1:4, lwd=0.75, lty=3, col=rgb(0.5, 0.5, 0.5, 0.6))
segments(x0=1:4, y0=g.o.b.l.c.2.age.pred.l95, y1=g.o.b.l.c.2.age.pred.u95, lwd=1.5, 
         col=magma(1, begin=0.3, end=0.3))
points(x=1:4, y=g.o.b.l.c.2.age.pred.mean, pch=16, col=magma(1, begin=0.3, end=0.3), cex=1.2)
text(x=4.2, y=0.85, labels="AmericasBarometer", pos=2, cex=1, adj=1)
box()
dev.off()

# year plot

year.pred = c(2006, 2008, 2010, 2012, 2014, 2017, 2019)
pred.dat.year = data.frame(birthyr=1964, republican=0, democrat=0, degree=0, female=1, white=1, income_std=0, 
                           south=0, year=year.pred, agecat4=3)
g.o.b.l.c.2.year.pred.eprob = posterior_epred(g.o.b.l.c.2.par, newdata=pred.dat.year, ndraws=1000)
g.o.b.l.c.2.year.pred.eprob2 = g.o.b.l.c.2.year.pred.eprob[,,5] + g.o.b.l.c.2.year.pred.eprob[,,6] + g.o.b.l.c.2.year.pred.eprob[,,7]
g.o.b.l.c.2.year.pred.mean = apply(g.o.b.l.c.2.year.pred.eprob2, 2, mean)
g.o.b.l.c.2.year.pred.l95 = apply(g.o.b.l.c.2.year.pred.eprob2, 2, quantile, probs=0.025)
g.o.b.l.c.2.year.pred.u95 = apply(g.o.b.l.c.2.year.pred.eprob2, 2, quantile, probs=0.975)
year_labs = c("2006", "2008", "2010", "2012", "2014", "2017", "2019")

pdf("fig4c.pdf", height=4, width=5.5)
par(mfrow=c(1,1), mar=c(3, 3, 2.2, 0.5), tcl=-0.2, las=0, cex=1)
plot(x=1:7, y=g.o.b.l.c.2.year.pred.mean, type="n", axes=FALSE, xlab="", ylab="", xlim=c(0.8, 7.2), 
     ylim=c(0.5, 1), yaxs="i")
axis(side=1, at=1:7, labels=year_labs, mgp=c(1, 0.1, 0), cex.axis=0.99)
axis(side=2, at=seq(0.5, 1, by=0.1), mgp=c(1, 0.3, 0), cex.axis=0.99)
mtext(side=1, line=1.3, "Year", cex=1.1)
mtext(side=2, line=1.4, "Probability of agreeing democracy is best (5-7)", cex=1)
mtext(side=3, line=0.5, "Support for democracy by year, controls", cex=1.2, adj=0.5)
abline(h=seq(0.5, 0.9, by=0.1), lwd=0.75, lty=3, col=rgb(0.5, 0.5, 0.5, 0.7))
abline(h=0, lwd=0.75, col=rgb(0.5, 0.5, 0.5, 0.6))
abline(v=1:7, lwd=0.75, lty=3, col=rgb(0.5, 0.5, 0.5, 0.6))
segments(x0=1:7, y0=g.o.b.l.c.2.year.pred.l95, y1=g.o.b.l.c.2.year.pred.u95, lwd=1.5, 
         col=magma(1, begin=0.3, end=0.3))
points(x=1:7, y=g.o.b.l.c.2.year.pred.mean, pch=16, col=magma(1, begin=0.3, end=0.3), cex=1.2)
text(x=7.2, y=0.98, labels="AmericasBarometer", pos=2, cex=1, adj=1)
box()
dev.off()

# cohort plot

g.o.b.l.c.2.coh.eff = conditional_effects(g.o.b.l.c.2.par, "birthyr", categorical=TRUE)

pred.dat = data.frame(birthyr=1913:2000, republican=0, democrat=0, degree=0, female=1, white=1, income_std=0, 
                      south=0, year=2012, agecat4=2)
g.o.b.l.c.2.pred.eprob = posterior_epred(g.o.b.l.c.2.par, newdata=pred.dat, ndraws=1000)
g.o.b.l.c.2.pred.eprob2 = g.o.b.l.c.2.pred.eprob[,,5] + g.o.b.l.c.2.pred.eprob[,,6] + g.o.b.l.c.2.pred.eprob[,,7]
g.o.b.l.c.2.pred.mean = apply(g.o.b.l.c.2.pred.eprob2, 2, mean)
g.o.b.l.c.2.pred.l95 = apply(g.o.b.l.c.2.pred.eprob2, 2, quantile, probs=0.025)
g.o.b.l.c.2.pred.u95 = apply(g.o.b.l.c.2.pred.eprob2, 2, quantile, probs=0.975)

pdf("fig2c.pdf", height=4, width=5.5)
par(mfrow=c(1,1), mar=c(3, 3, 2.2, 0.5), tcl=-0.2, las=0, cex=1)
plot(x=1913:2000, y=g.o.b.l.c.2.pred.mean, type="n", axes=FALSE, xlab="", ylab="", xaxs="i", yaxs="i",
     ylim=c(0.3, 1))
axis(side=1, at=seq(1910, 1990, by=10), mgp=c(1, 0.1, 0), cex.axis=0.99)
axis(side=2, at=seq(0.4, 1, by=0.2), mgp=c(1, 0.3, 0), cex.axis=0.99)
mtext(side=1, line=1.3, "Year of birth", cex=1.1)
mtext(side=2, line=1.4, "Probability of agreeing democracy is best (5-7)", cex=1)
mtext(side=3, line=0.5, "Support for democracy by birth year, controls", cex=1.2, adj=0.5)
abline(h=seq(0.4, 1, by=0.2), lwd=0.75, lty=3, col=rgb(0.5, 0.5, 0.5, 0.7))
abline(h=0, lwd=0.75, col=rgb(0.5, 0.5, 0.5, 0.6))
abline(v=seq(1920, 1990, by=10), lwd=0.75, lty=3, col=rgb(0.5, 0.5, 0.5, 0.6))
lines(x=1913:2000, y=g.o.b.l.c.2.pred.mean, lwd=2, col=magma(1, begin=0.3, end=0.3))
polygon(x=c(1913:2000, 2000:1913), y=c(g.o.b.l.c.2.pred.l95, g.o.b.l.c.2.pred.u95[88:1]), 
        col=magma(1, begin=0.3, end=0.3, alpha=0.5), border=NA, density=NA)
text(x=2000, y=0.96, labels="AmericasBarometer", pos=2, cex=1)
box()
dev.off()


## Bayesian ordinal GAMs using BRMS - democracy is best / LAPOP - no controls - in paper

lapop2 = lapop
lapop2$year = as.factor(lapop$year)
lapop2$agecat4 = as.factor(lapop$agecat4)
lapop2$dem_best = as.ordered(lapop$dem_best)

g.o.b.l.nc.2.ser <- brm(dem_best | weights(wt) ~ s(birthyr) + year + agecat4, data=lapop2, prior=g.o.b.l.prior, 
                        family=cumulative("logit"), control=list(adapt_delta=0.95), chains=6, iter=12, 
                        backend="cmdstanr")

g.o.b.l.nc.2.par <- update(
  g.o.b.l.nc.2.ser, chains=3, cores=3, iter=2000, backend="cmdstanr", threads=threading(4),
  control=c(adapt_delta=0.95, max_treedepth=12), stan_model_args=list(stanc_options = list("O1"))
)

summary(g.o.b.l.nc.2.par)
pp_check(g.o.b.l.nc.2.par)
pp_check(g.o.b.l.nc.2.par, type="ecdf_overlay")

pdf("fig_s7a.pdf", width=6, height=6)
plot(g.o.b.l.nc.2.par, newpage=FALSE, 
     variable=c("b_Intercept[1]", "b_year2008", "b_agecat42", "bs_sbirthyr_1", "sds_sbirthyr_1"))
dev.off()

sum.g.o.b.l.nc.2.par <- capture.output(summary(g.o.b.l.nc.2.par))
cat(title="Lapop GAM results", sum.g.o.b.l.nc.2.par, file="table2_1.txt", sep="\n", append=TRUE)

# age plot

pred.dat.age = data.frame(birthyr=1964, year=2012, agecat4=1:4)
g.o.b.l.nc.2.age.pred.eprob = posterior_epred(g.o.b.l.nc.2.par, newdata=pred.dat.age, ndraws=1000)
g.o.b.l.nc.2.age.pred.eprob2 = g.o.b.l.nc.2.age.pred.eprob[,,5] + g.o.b.l.nc.2.age.pred.eprob[,,6] + g.o.b.l.nc.2.age.pred.eprob[,,7]
g.o.b.l.nc.2.age.pred.mean = apply(g.o.b.l.nc.2.age.pred.eprob2, 2, mean)
g.o.b.l.nc.2.age.pred.l95 = apply(g.o.b.l.nc.2.age.pred.eprob2, 2, quantile, probs=0.025)
g.o.b.l.nc.2.age.pred.u95 = apply(g.o.b.l.nc.2.age.pred.eprob2, 2, quantile, probs=0.975)

pdf("fig3a.pdf", height=4, width=5.5)
par(mfrow=c(1,1), mar=c(3, 3, 2.2, 0.5), tcl=-0.2, las=0, cex=1)
plot(x=1:4, y=g.o.b.l.nc.2.age.pred.mean, type="n", axes=FALSE, xlab="", ylab="", xlim=c(0.8, 4.2), 
     ylim=c(0.6, 0.85))
axis(side=1, at=1:4, labels=c("18-29", "30-44", "45-59", "60+"), mgp=c(1, 0.1, 0), cex.axis=0.99)
axis(side=2, at=seq(0.6, 0.9, by=0.1), mgp=c(1, 0.3, 0), cex.axis=0.99)
mtext(side=1, line=1.3, "Age group", cex=1.1)
mtext(side=2, line=1.4, "Probability of agreeing democracy is best (5-7)", cex=1)
mtext(side=3, line=0.5, "Support for democracy by age group, no controls", cex=1.2, adj=0.5)
abline(h=seq(0.7, 0.9, by=0.05), lwd=0.75, lty=3, col=rgb(0.5, 0.5, 0.5, 0.7))
abline(h=0, lwd=0.75, col=rgb(0.5, 0.5, 0.5, 0.6))
abline(v=1:4, lwd=0.75, lty=3, col=rgb(0.5, 0.5, 0.5, 0.6))
segments(x0=1:4, y0=g.o.b.l.nc.2.age.pred.l95, y1=g.o.b.l.nc.2.age.pred.u95, lwd=1.5, 
         col=magma(1, begin=0.3, end=0.3))
points(x=1:4, y=g.o.b.l.nc.2.age.pred.mean, pch=16, col=magma(1, begin=0.3, end=0.3), cex=1.2)
text(x=4.2, y=0.85, labels="AmericasBarometer", pos=2, cex=1, adj=1)
box()
dev.off()

# year plot

year.pred = c(2006, 2008, 2010, 2012, 2014, 2017, 2019)
pred.dat.year = data.frame(birthyr=1964, year=year.pred, agecat4=3)
g.o.b.l.nc.2.year.pred.eprob = posterior_epred(g.o.b.l.nc.2.par, newdata=pred.dat.year, ndraws=1000)
g.o.b.l.nc.2.year.pred.eprob2 = g.o.b.l.nc.2.year.pred.eprob[,,5] + g.o.b.l.nc.2.year.pred.eprob[,,6] + g.o.b.l.nc.2.year.pred.eprob[,,7]
g.o.b.l.nc.2.year.pred.mean = apply(g.o.b.l.nc.2.year.pred.eprob2, 2, mean)
g.o.b.l.nc.2.year.pred.l95 = apply(g.o.b.l.nc.2.year.pred.eprob2, 2, quantile, probs=0.025)
g.o.b.l.nc.2.year.pred.u95 = apply(g.o.b.l.nc.2.year.pred.eprob2, 2, quantile, probs=0.975)
year_labs = c("2006", "2008", "2010", "2012", "2014", "2017", "2019")

pdf("fig4a.pdf", height=4, width=5.5)
par(mfrow=c(1,1), mar=c(3, 3, 2.2, 0.5), tcl=-0.2, las=0, cex=1)
plot(x=1:7, y=g.o.b.l.nc.2.year.pred.mean, type="n", axes=FALSE, xlab="", ylab="", xlim=c(0.8, 7.2), 
     ylim=c(0.5, 1), yaxs="i")
axis(side=1, at=1:7, labels=year_labs, mgp=c(1, 0.1, 0), cex.axis=0.99)
axis(side=2, at=seq(0.5, 1, by=0.1), mgp=c(1, 0.3, 0), cex.axis=0.99)
mtext(side=1, line=1.3, "Year", cex=1.1)
mtext(side=2, line=1.4, "Probability of agreeing democracy is best (5-7)", cex=1)
mtext(side=3, line=0.5, "Support for democracy by year, no controls", cex=1.2, adj=0.5)
abline(h=seq(0.5, 1, by=0.1), lwd=0.75, lty=3, col=rgb(0.5, 0.5, 0.5, 0.7))
abline(h=0, lwd=0.75, col=rgb(0.5, 0.5, 0.5, 0.6))
abline(v=1:7, lwd=0.75, lty=3, col=rgb(0.5, 0.5, 0.5, 0.6))
segments(x0=1:7, y0=g.o.b.l.nc.2.year.pred.l95, y1=g.o.b.l.nc.2.year.pred.u95, lwd=1.5, 
         col=magma(1, begin=0.3, end=0.3))
points(x=1:7, y=g.o.b.l.nc.2.year.pred.mean, pch=16, col=magma(1, begin=0.3, end=0.3), cex=1.2)
text(x=7.2, y=0.98, labels="AmericasBarometer", pos=2, cex=1, adj=1)
box()
dev.off()

# cohort plot

g.o.b.l.nc.2.coh.eff = conditional_effects(g.o.b.l.nc.2.par, "birthyr", categorical=TRUE)

pred.dat = data.frame(birthyr=1913:2000, year=2012, agecat4=2)
g.o.b.l.nc.2.pred.eprob = posterior_epred(g.o.b.l.nc.2.par, newdata=pred.dat, ndraws=1000)
g.o.b.l.nc.2.pred.eprob2 = g.o.b.l.nc.2.pred.eprob[,,5] + g.o.b.l.nc.2.pred.eprob[,,6] + g.o.b.l.nc.2.pred.eprob[,,7]
g.o.b.l.nc.2.pred.mean = apply(g.o.b.l.nc.2.pred.eprob2, 2, mean)
g.o.b.l.nc.2.pred.l95 = apply(g.o.b.l.nc.2.pred.eprob2, 2, quantile, probs=0.025)
g.o.b.l.nc.2.pred.u95 = apply(g.o.b.l.nc.2.pred.eprob2, 2, quantile, probs=0.975)

pdf("fig2a.pdf", height=4, width=5.5)
par(mfrow=c(1,1), mar=c(3, 3, 2.2, 0.5), tcl=-0.2, las=0, cex=1)
plot(x=1913:2000, y=g.o.b.l.nc.2.pred.mean, type="n", axes=FALSE, xlab="", ylab="", xaxs="i", yaxs="i",
     ylim=c(0.3, 1))
axis(side=1, at=seq(1910, 1990, by=10), mgp=c(1, 0.1, 0), cex.axis=0.99)
axis(side=2, at=seq(0.4, 1, by=0.2), mgp=c(1, 0.3, 0), cex.axis=0.99)
mtext(side=1, line=1.3, "Year of birth", cex=1.1)
mtext(side=2, line=1.4, "Probability of agreeing democracy is best (5-7)", cex=1)
mtext(side=3, line=0.5, "Support for democracy by birth year, no controls", cex=1.2, adj=0.5)
abline(h=seq(0.4, 1, by=0.1), lwd=0.75, lty=3, col=rgb(0.5, 0.5, 0.5, 0.7))
abline(h=0, lwd=0.75, col=rgb(0.5, 0.5, 0.5, 0.6))
abline(v=seq(1920, 1990, by=10), lwd=0.75, lty=3, col=rgb(0.5, 0.5, 0.5, 0.6))
lines(x=1913:2000, y=g.o.b.l.nc.2.pred.mean, lwd=2, col=magma(1, begin=0.3, end=0.3))
polygon(x=c(1913:2000, 2000:1913), y=c(g.o.b.l.nc.2.pred.l95, g.o.b.l.nc.2.pred.u95[88:1]), 
        col=magma(1, begin=0.3, end=0.3, alpha=0.5), border=NA, density=NA)
text(x=2000, y=0.96, labels="AmericasBarometer", pos=2, cex=1)
box()
dev.off()



### WVS

## Bayesian HAPC - WVS - w / controls - in SI

h.b.w.c.1.prior = c(set_prior("normal(0, 1)", class="b"),
                    set_prior("normal(0, 2)", class="Intercept"),
                    set_prior("normal(0, 1)", class="sd"),
                    set_prior("normal(0, 2)", class="sigma"))

h.b.w.c.1.ser <- brm(sup_dem | weights(wt) ~ republican + democrat + female + white + income_std + degree + south
                     + (1 | year) + (1 | agecat4) + (1 | decbirth), data=wvs, 
                     prior=h.b.w.c.1.prior, control=list(adapt_delta=0.99), chains=9, iter=12, backend="cmdstanr")

h.b.w.c.1.par <- update(
  h.b.w.c.1.ser, chains=3, cores=3, iter=2000, backend="cmdstanr", threads=threading(4),
  control=c(adapt_delta=0.99, max_treedepth=12), stan_model_args=list(stanc_options = list("O1"))
)

summary(h.b.w.c.1.par)
pp_check(h.b.w.c.1.par)

pdf("fig_s10b.pdf", width=6, height=6)
plot(h.b.w.c.1.par, newpage=FALSE, 
     variable=c("b_Intercept", "sd_agecat4__Intercept", "sd_decbirth__Intercept", "sd_year__Intercept", "sigma"))
dev.off()

sum.h.b.w.c.1.par <- capture.output(summary(h.b.w.c.1.par))
cat(title="WVS HAPCM results", sum.h.b.w.c.1.par, file="table_s3_2.txt", sep="\n", append=TRUE)

# cohort plots

gen_pe = ranef(h.b.w.c.1.par)$decbirth[,,1][,1]
gen_up = ranef(h.b.w.c.1.par)$decbirth[,,1][,4]
gen_lo = ranef(h.b.w.c.1.par)$decbirth[,,1][,3]
gen_labs = c("1900s", "1910s", "1920s", "1930s", "1940s", "1950s", "1960s", "1970s", "1980s", "1990s")

pdf("fig_s3d.pdf", height=4, width=5.5)
plot.col = viridis(n=1, begin=0.4, end=0.4)
plot.alp = viridis(n=1, begin=0.4, end=0.4, alpha=0.5)
par(mfrow=c(1,1), mar=c(3, 3, 2.2, 0.5), tcl=-0.2, las=0)
plot(x=1:10, y=gen_pe, pch=16, xlim=c(1, 10), ylim=c(-0.4, 0.4), type="n",
     xlab="", ylab="", main="", axes=FALSE)
abline(h=seq(-0.4, 0.4, by=0.2), col=rgb(0.5, 0.5, 0.5, 0.8), lwd=0.75, lty=3)
abline(v=1:10, col=rgb(0.5, 0.5, 0.5, 0.8), lwd=0.75, lty=3)
lines(x=1:10, y=gen_pe, lwd=2, col=plot.col)
polygon(x=c(1:10, 10:1), y=c(gen_up, gen_lo[10:1]), col=plot.alp, border=NA, density=NA)
axis(1, at=1:10, labels=gen_labs, cex.axis=0.99, las=1, mgp=c(1, 0.15, 0))
axis(2, at=seq(-0.4, 0.4, by=0.2), cex.axis=0.99, las=1, mgp=c(1, 0.35, 0))
mtext(side=1, line=1.3, "Decade of birth", cex=1.1)
mtext(side=2, line=1.5, "Support democracy & rejection autocracy, std scale", cex=1)
mtext("Support for democracy by birth decade, controls", side=3, line=0.5, adj=0.5, cex=1.2)
text(x=10, y=0.4, labels="WVS", pos=2, cex=1)
box()
dev.off()

# age plots

age_pe = ranef(h.b.w.c.1.par)$agecat4[,,1][,1]
age_up = ranef(h.b.w.c.1.par)$agecat4[,,1][,4]
age_lo = ranef(h.b.w.c.1.par)$agecat4[,,1][,3]
age_labs = c("18-29", "30-44", "45-59", "60+")

pdf("fig_s4d.pdf", height=4, width=5.5)
plot.col = viridis(n=1, begin=0.4, end=0.4)
plot.alp = viridis(n=1, begin=0.4, end=0.4, alpha=0.5)
par(mfrow=c(1,1), mar=c(3, 3, 2.2, 0.5), tcl=-0.2, las=0)
plot(x=1:4, y=age_pe, pch=16, xlim=c(0.8, 4.2), ylim=c(-0.6, 0.7), type="n",
     xlab="", ylab="", main="", axes=FALSE)
abline(h=seq(-1, 1.5, by=0.5), col=rgb(0.5, 0.5, 0.5, 0.9), lwd=0.75, lty=3)
abline(v=1:10, col=rgb(0.5, 0.5, 0.5, 0.9), lwd=0.75, lty=3)
points(x=1:4, y=age_pe, pch=16, cex=1.2, col=plot.col)
segments(x0=1:4, y0=age_lo, y1=age_up, lwd=1.5, col=plot.col)
axis(1, at=1:4, labels=age_labs, cex.axis=0.99, las=1, mgp=c(1, 0.15, 0))
axis(2, at=seq(-1, 1.5, by=0.5), cex.axis=0.99, las=1, mgp=c(1, 0.35, 0))
text(x=4.2, y=0.7, labels="WVS", pos=2, cex=1, adj=1)
mtext(side=1, line=1.3, "Age group", cex=1.1)
mtext(side=2, line=1.5, "Support democracy & rejection autocracy, std scale", cex=1)
mtext("Support for democracy by age, controls", side=3, line=0.5, adj=0.5, cex=1.2)
box()
dev.off()

# year plots

year_pe = ranef(h.b.w.c.1.par)$year[,,1][,1]
year_up = ranef(h.b.w.c.1.par)$year[,,1][,4]
year_lo = ranef(h.b.w.c.1.par)$year[,,1][,3]
year_labs = c("1995", "1999", "2006", "2011", "2017")

pdf("fig_s5d.pdf", height=4, width=5.5)
plot.col = viridis(n=1, begin=0.4, end=0.4)
plot.alp = viridis(n=1, begin=0.4, end=0.4, alpha=0.5)
par(mfrow=c(1,1), mar=c(3, 3, 2.2, 0.5), tcl=-0.2, las=0)
plot(x=1:5, y=year_pe, pch=16, xlim=c(0.8, 5.2), ylim=c(-0.6, 0.7), type="n",
     xlab="", ylab="", main="", axes=FALSE)
abline(h=seq(-1, 1.5, by=0.5), col=rgb(0.5, 0.5, 0.5, 0.9), lwd=0.75, lty=3)
abline(v=1:10, col=rgb(0.5, 0.5, 0.5, 0.9), lwd=0.75, lty=3)
points(x=1:5, y=year_pe, pch=16, cex=1.2, col=plot.col)
segments(x0=1:5, y0=year_lo, y1=year_up, lwd=1.5, col=plot.col)
axis(1, at=1:5, labels=year_labs, cex.axis=0.99, las=1, mgp=c(1, 0.15, 0))
axis(2, at=seq(-1, 1.5, by=0.5), cex.axis=0.99, las=1, mgp=c(1, 0.35, 0))
text(x=5.2, y=0.7, labels="WVS", pos=2, cex=1, adj=1)
mtext(side=1, line=1.3, "Year", cex=1.1)
mtext(side=2, line=1.6, "Support democracy & rejection autocracy, std scale", cex=1)
mtext("Support for democracy by year, controls", side=3, line=0.5, adj=0.5, cex=1.2)
box()
dev.off()


## Bayesian HAPC - WVS - no controls - in SI

h.b.w.nc.1.prior = c(set_prior("normal(0, 2)", class="Intercept"),
                     set_prior("normal(0, 1)", class="sd"),
                     set_prior("normal(0, 2)", class="sigma"))

h.b.w.nc.1.ser <- brm(sup_dem | weights(wt) ~ (1 | year) + (1 | agecat4) + (1 | decbirth), 
                      data=wvs, prior=h.b.w.nc.1.prior, control=list(adapt_delta=0.99), 
                      chains=9, iter=12, backend="cmdstanr")

h.b.w.nc.1.par <- update(
  h.b.w.nc.1.ser, chains=3, cores=3, iter=2000, backend="cmdstanr", threads=threading(4),
  control=c(adapt_delta=0.99, max_treedepth=12), stan_model_args=list(stanc_options = list("O1"))
)

summary(h.b.w.nc.1.par)
pp_check(h.b.w.nc.1.par)

pdf("fig_s10a.pdf", width=6, height=6)
plot(h.b.w.nc.1.par, newpage=FALSE, 
     variable=c("b_Intercept", "sd_agecat4__Intercept", "sd_decbirth__Intercept", "sd_year__Intercept", "sigma"))
dev.off()

sum.h.b.w.nc.1.par <- capture.output(summary(h.b.w.nc.1.par))
cat(title="WVS HAPCM results", sum.h.b.w.nc.1.par, file="table_s2_1.txt", sep="\n", append=TRUE)

# cohort plots

gen_pe = ranef(h.b.w.nc.1.par)$decbirth[,,1][,1]
gen_up = ranef(h.b.w.nc.1.par)$decbirth[,,1][,4]
gen_lo = ranef(h.b.w.nc.1.par)$decbirth[,,1][,3]
gen_labs = c("1900s", "1910s", "1920s", "1930s", "1940s", "1950s", "1960s", "1970s", "1980s", "1990s")

pdf("fig_s3b.pdf", height=4, width=5.5)
plot.col = viridis(n=1, begin=0.4, end=0.4)
plot.alp = viridis(n=1, begin=0.4, end=0.4, alpha=0.5)
par(mfrow=c(1,1), mar=c(3, 3, 2.2, 0.5), tcl=-0.2, las=0)
plot(x=1:10, y=gen_pe, pch=16, xlim=c(1, 10), ylim=c(-0.4, 0.4), type="n",
     xlab="", ylab="", main="", axes=FALSE)
abline(h=seq(-0.4, 0.4, by=0.2), col=rgb(0.5, 0.5, 0.5, 0.8), lwd=0.75, lty=3)
abline(v=1:10, col=rgb(0.5, 0.5, 0.5, 0.8), lwd=0.75, lty=3)
lines(x=1:10, y=gen_pe, lwd=2, col=plot.col)
polygon(x=c(1:10, 10:1), y=c(gen_up, gen_lo[10:1]), col=plot.alp, border=NA, density=NA)
axis(1, at=1:10, labels=gen_labs, cex.axis=0.99, las=1, mgp=c(1, 0.15, 0))
axis(2, at=seq(-0.4, 0.4, by=0.2), cex.axis=0.99, las=1, mgp=c(1, 0.35, 0))
mtext(side=1, line=1.3, "Decade of birth", cex=1.1)
mtext(side=2, line=1.5, "Support democracy & rejection autocracy, std scale", cex=1)
mtext("Support for democracy by birth decade, no controls", side=3, line=0.5, adj=0.5, cex=1.2)
text(x=10, y=0.4, labels="WVS", pos=2, cex=1)
box()
dev.off()

# age plots

age_pe = ranef(h.b.w.nc.1.par)$agecat4[,,1][,1]
age_up = ranef(h.b.w.nc.1.par)$agecat4[,,1][,4]
age_lo = ranef(h.b.w.nc.1.par)$agecat4[,,1][,3]
age_labs = c("18-29", "30-44", "45-59", "60+")

pdf("fig_s4b.pdf", height=4, width=5.5)
plot.col = viridis(n=1, begin=0.4, end=0.4)
plot.alp = viridis(n=1, begin=0.4, end=0.4, alpha=0.5)
par(mfrow=c(1,1), mar=c(3, 3, 2.2, 0.5), tcl=-0.2, las=0)
plot(x=1:4, y=age_pe, pch=16, xlim=c(0.8, 4.2), ylim=c(-0.6, 0.7), type="n",
     xlab="", ylab="", main="", axes=FALSE)
abline(h=seq(-1, 1.5, by=0.5), col=rgb(0.5, 0.5, 0.5, 0.9), lwd=0.75, lty=3)
abline(v=1:10, col=rgb(0.5, 0.5, 0.5, 0.9), lwd=0.75, lty=3)
points(x=1:4, y=age_pe, pch=16, cex=1.2, col=plot.col)
segments(x0=1:4, y0=age_lo, y1=age_up, lwd=1.5, col=plot.col)
axis(1, at=1:4, labels=age_labs, cex.axis=0.99, las=1, mgp=c(1, 0.15, 0))
axis(2, at=seq(-1, 1.5, by=0.5), cex.axis=0.99, las=1, mgp=c(1, 0.35, 0))
text(x=4.2, y=0.7, labels="WVS", pos=2, cex=1, adj=1)
mtext(side=1, line=1.3, "Age group", cex=1.1)
mtext(side=2, line=1.5, "Support democracy & rejection autocracy, std scale", cex=1)
mtext("Support for democracy by age, no controls", side=3, line=0.5, adj=0.5, cex=1.2)
box()
dev.off()

# year plots

year_pe = ranef(h.b.w.nc.1.par)$year[,,1][,1]
year_up = ranef(h.b.w.nc.1.par)$year[,,1][,4]
year_lo = ranef(h.b.w.nc.1.par)$year[,,1][,3]
year_labs = c("1995", "1999", "2006", "2011", "2017")

pdf("fig_s5b.pdf", height=4, width=5.5)
plot.col = viridis(n=1, begin=0.4, end=0.4)
plot.alp = viridis(n=1, begin=0.4, end=0.4, alpha=0.5)
par(mfrow=c(1,1), mar=c(3, 3, 2.2, 0.5), tcl=-0.2, las=0)
plot(x=1:5, y=year_pe, pch=16, xlim=c(0.8, 5.2), ylim=c(-0.6, 0.7), type="n",
     xlab="", ylab="", main="", axes=FALSE)
abline(h=seq(-1, 1.5, by=0.5), col=rgb(0.5, 0.5, 0.5, 0.9), lwd=0.75, lty=3)
abline(v=1:10, col=rgb(0.5, 0.5, 0.5, 0.9), lwd=0.75, lty=3)
points(x=1:5, y=year_pe, pch=16, cex=1.2, col=plot.col)
segments(x0=1:5, y0=year_lo, y1=year_up, lwd=1.5, col=plot.col)
axis(1, at=1:5, labels=year_labs, cex.axis=0.99, las=1, mgp=c(1, 0.15, 0))
axis(2, at=seq(-1, 1.5, by=0.5), cex.axis=0.99, las=1, mgp=c(1, 0.35, 0))
text(x=5.2, y=0.7, labels="WVS", pos=2, cex=1, adj=1)
mtext(side=1, line=1.3, "Year", cex=1.1)
mtext(side=2, line=1.6, "Support democracy & rejection autocracy, std scale", cex=1)
mtext("Support for democracy by year, no controls", side=3, line=0.5, adj=0.5, cex=1.2)
box()
dev.off()


## Bayesian GAMMs - WVS - 3-item Rej Autoc - 4 age groups

g.b.w.prior = c(set_prior("normal(0, 1)", class="b"),
                set_prior("normal(0, 1)", class="Intercept"),
                set_prior("normal(0, 2)", class="sigma"),
                set_prior("normal(0, 2)", class="sds"))

wvs2 = wvs
wvs2$year = as.factor(wvs2$year)
wvs2$agecat4 = as.factor(wvs2$agecat4)

g.b.w.c.3.ser = brm(rej_aut | weights(wt) ~ s(birthyr) + republican + democrat + degree + female + white + south
                    + income_std + year + agecat4, data=wvs2, prior=g.b.w.prior,
                    control=list(adapt_delta=0.99), chains=6, iter=12, backend="cmdstanr")

g.b.w.c.3.par <- update(
  g.b.w.c.3.ser, chains=3, cores=3, iter=2000, backend="cmdstanr", threads=threading(4),
   control=c(adapt_delta=0.99, max_treedepth=12), stan_model_args=list(stanc_options = list("O1"))
)

summary(g.b.w.c.3.par)
pp_check(g.b.w.c.3.par)
pp_check(g.b.w.c.3.par, type="ecdf_overlay")

plot(g.b.w.c.3.par, newpage=FALSE, 
     variable=c("b_year1999", "b_agecat42", "bs_sbirthyr_1", "sds_sbirthyr_1", "sigma"))

# generational effects

g.b.w.c.3.coh.eff = conditional_effects(g.b.w.c.3.par, "birthyr")

pdf("fig_s6a.pdf", height=4, width=5.5)
par(mfrow=c(1,1), mar=c(3, 3, 2.2, 0.5), tcl=-0.2, las=0, cex=1)
plot(x=g.b.w.c.3.coh.eff[[1]]$birthyr, y=g.b.w.c.3.coh.eff[[1]]$estimate, type="n", axes=FALSE,
     xlab="", ylab="", ylim=c(-0.3, 0.8), xaxs="i")
axis(side=1, at=seq(1900, 2000, by=10), mgp=c(1, 0.1, 0), cex.axis=0.99)
axis(side=2, at=seq(-0.2, 0.8, by=0.2), mgp=c(1, 0.3, 0), cex.axis=0.99)
mtext(side=1, line=1.3, "Year of birth", cex=1.1)
mtext(side=2, line=1.4, "Reject autocracy (3-item), std scale", cex=1)
mtext(side=3, line=0.5, "Reject autocracy by birth year, GAM", cex=1.2, adj=0.5)
abline(h=seq(-0.2, 0.8, by=0.2), lwd=0.75, lty=3, col=rgb(0.5, 0.5, 0.5, 0.7))
abline(v=seq(1910, 1990, by=10), lwd=0.75, lty=3, col=rgb(0.5, 0.5, 0.5, 0.7))
lines(x=g.b.w.c.3.coh.eff[[1]]$birthyr, y=g.b.w.c.3.coh.eff[[1]]$estimate, lwd=2,
      col=magma(1, begin=0.3, end=0.3))
polygon(x=c(g.b.w.c.3.coh.eff[[1]]$birthyr, g.b.w.c.3.coh.eff[[1]]$birthyr[100:1]),
        y=c(g.b.w.c.3.coh.eff[[1]]$lower, g.b.w.c.3.coh.eff[[1]]$upper[100:1]),
        col=magma(1, begin=0.3, end=0.3, alpha=0.5), border=NA, density=NA)
text(x=2000, y=0.8, labels="WVS", pos=2, cex=1)
box()
dev.off()


## Bayesian GAMMs - WVS - 4-item Sup Dem - 4 age groups - controls

g.b.w.prior = c(set_prior("normal(0, 1)", class="b"),
                set_prior("normal(0, 1)", class="Intercept"),
                set_prior("normal(0, 2)", class="sigma"),
                set_prior("normal(0, 2)", class="sds"))

wvs2 = wvs
wvs2$year = as.factor(wvs2$year)
wvs2$agecat4 = as.factor(wvs2$agecat4)

g.b.w.c.4.ser = brm(sup_dem_2 | weights(wt) ~ s(birthyr) + republican + democrat + degree + female + white + south
                    + income_std + year + agecat4, data=wvs2, prior=g.b.w.prior,
                    control=list(adapt_delta=0.99), chains=6, iter=12, backend="cmdstanr")

g.b.w.c.4.par <- update(
  g.b.w.c.4.ser, chains=3, cores=3, iter=2000, backend="cmdstanr", threads=threading(4),
  control=c(adapt_delta=0.99, max_treedepth=12), stan_model_args=list(stanc_options = list("O1"))
)

summary(g.b.w.c.4.par)
pp_check(g.b.w.c.4.par)
pp_check(g.b.w.c.4.par, type="ecdf_overlay")

plot(g.b.w.c.4.par, newpage=FALSE, 
     variable=c("b_year1999", "b_agecat42", "bs_sbirthyr_1", "sds_sbirthyr_1", "sigma"))

# generational effects

g.b.w.c.4.coh.eff = conditional_effects(g.b.w.c.4.par, "birthyr")

pdf("fig_s6b.pdf", height=4, width=5.5)
par(mfrow=c(1,1), mar=c(3, 3, 2.2, 0.5), tcl=-0.2, las=0, cex=1)
plot(x=g.b.w.c.4.coh.eff[[1]]$birthyr, y=g.b.w.c.4.coh.eff[[1]]$estimate, type="n", axes=FALSE,
     xlab="", ylab="", ylim=c(-0.3, 0.8), xaxs="i")
axis(side=1, at=seq(1900, 2000, by=10), mgp=c(1, 0.1, 0), cex.axis=0.99)
axis(side=2, at=seq(-0.2, 0.8, by=0.2), mgp=c(1, 0.3, 0), cex.axis=0.99)
mtext(side=1, line=1.3, "Year of birth", cex=1.1)
mtext(side=2, line=1.4, "Support democracy (4-item), std scale", cex=1)
mtext(side=3, line=0.5, "Support democracy by birth year, GAM", cex=1.2, adj=0.5)
abline(h=seq(-0.2, 0.8, by=0.2), lwd=0.75, lty=3, col=rgb(0.5, 0.5, 0.5, 0.7))
abline(v=seq(1910, 1990, by=10), lwd=0.75, lty=3, col=rgb(0.5, 0.5, 0.5, 0.7))
lines(x=g.b.w.c.4.coh.eff[[1]]$birthyr, y=g.b.w.c.4.coh.eff[[1]]$estimate, lwd=2,
      col=magma(1, begin=0.3, end=0.3))
polygon(x=c(g.b.w.c.4.coh.eff[[1]]$birthyr, g.b.w.c.4.coh.eff[[1]]$birthyr[100:1]),
        y=c(g.b.w.c.4.coh.eff[[1]]$lower, g.b.w.c.4.coh.eff[[1]]$upper[100:1]),
        col=magma(1, begin=0.3, end=0.3, alpha=0.5), border=NA, density=NA)
text(x=2000, y=0.8, labels="WVS", pos=2, cex=1)
box()
dev.off()



## Bayesian GAMs - support democracy & reject autocracy - no controls - in paper

g.b.w.2.prior = c(set_prior("normal(0, 1)", class="b"),
                  set_prior("normal(0, 1)", class="Intercept"),
                  set_prior("normal(0, 2)", class="sigma"),
                  set_prior("normal(0, 2)", class="sds"))

wvs2 = wvs
wvs2$year = as.factor(wvs2$year)
wvs2$agecat4 = as.factor(wvs2$agecat4)
g.b.w.nc.2.ser = brm(sup_dem | weights(wt) ~ s(birthyr) + year + agecat4, data=wvs2, prior=g.b.w.2.prior, 
                     control=list(adapt_delta=0.99), chains=6, iter=12, backend="cmdstanr")

g.b.w.nc.2.par <- update(
  g.b.w.nc.2.ser, chains=3, cores=3, iter=2000, backend="cmdstanr", threads=threading(4),
  control=c(adapt_delta=0.99, max_treedepth=12), stan_model_args=list(stanc_options = list("O1"))
)

summary(g.b.w.nc.2.par)
pp_check(g.b.w.nc.2.par)
pp_check(g.b.w.nc.2.par, type="ecdf_overlay")

pdf("fig_s8a.pdf", width=6, height=6)
plot(g.b.w.nc.2.par, newpage=FALSE, 
     variable=c("b_year1999", "b_agecat42", "bs_sbirthyr_1", "sds_sbirthyr_1", "sigma"))
dev.off()

sum.g.b.w.nc.2.par <- capture.output(summary(g.b.w.nc.2.par))
cat(title="WVS GAM results", sum.g.b.w.nc.2.par, file="table3_1.txt", sep="\n", append=TRUE)

# age plot

pred.dat.age = data.frame(birthyr=1961, year=2006, agecat4=1:4)
g.b.w.nc.2.age.pred.eprob = posterior_epred(g.b.w.nc.2.par, newdata=pred.dat.age, ndraws=1000)
g.b.w.nc.2.age.pred.mean = apply(g.b.w.nc.2.age.pred.eprob, 2, mean)
g.b.w.nc.2.age.pred.l95 = apply(g.b.w.nc.2.age.pred.eprob, 2, quantile, probs=0.025)
g.b.w.nc.2.age.pred.u95 = apply(g.b.w.nc.2.age.pred.eprob, 2, quantile, probs=0.975)

pdf("fig3b.pdf", height=4, width=5.5)
par(mfrow=c(1,1), mar=c(3, 3, 2.2, 0.5), tcl=-0.2, las=0, cex=1)
plot(x=1:4, y=g.b.w.nc.2.age.pred.mean, type="n", axes=FALSE, xlab="", ylab="", 
     xlim=c(0.8, 4.2), ylim=c(-0.4, 0.2))
axis(side=1, at=1:4, labels=c("18-29", "30-44", "45-59", "60+"), mgp=c(1, 0.1, 0), cex.axis=0.99)
axis(side=2, at=seq(-0.4, 0.2, by=0.2), mgp=c(1, 0.3, 0), cex.axis=0.99)
mtext(side=1, line=1.3, "Age group", cex=1.1)
mtext(side=2, line=1.4, "Support democracy & reject autocracy, std scale", cex=1)
mtext(side=3, line=0.5, "Support for democracy by age group, no controls", cex=1.2, adj=0.5)
abline(h=seq(-0.2, 0.2, by=0.1), lwd=0.75, lty=3, col=rgb(0.5, 0.5, 0.5, 0.7))
abline(v=1:4, lwd=0.75, lty=3, col=rgb(0.5, 0.5, 0.5, 0.7))
points(x=1:4, y=g.b.w.nc.2.age.pred.mean, pch=16, 
       col=magma(1, begin=0.3, end=0.3), cex=1.2)
segments(x0=1:4, y0=g.b.w.nc.2.age.pred.l95, y1=g.b.w.nc.2.age.pred.u95, col=magma(1, begin=0.3, end=0.3), lwd=1.5)
text(x=4.2, y=0.2, labels="WVS", pos=2, cex=1, adj=1)
box()
dev.off()

# cohort plot

g.b.w.nc.2.coh.eff = conditional_effects(g.b.w.nc.2.par, "birthyr")

pdf("fig2b.pdf", height=4, width=5.5)
par(mfrow=c(1,1), mar=c(3, 3, 2.2, 0.5), tcl=-0.2, las=0, cex=1)
plot(x=g.b.w.nc.2.coh.eff[[1]]$birthyr, y=g.b.w.nc.2.coh.eff[[1]]$estimate, type="n", axes=FALSE, 
     xlab="", ylab="", ylim=c(-0.4, 0.9), xaxs="i")
axis(side=1, at=seq(1900, 2000, by=10), mgp=c(1, 0.1, 0), cex.axis=0.99)
axis(side=2, at=seq(-0.4, 0.8, by=0.4), mgp=c(1, 0.3, 0), cex.axis=0.99)
mtext(side=1, line=1.3, "Year of birth", cex=1.1)
mtext(side=2, line=1.4, "Support democracy & reject autocracy, std scale", cex=1)
mtext(side=3, line=0.5, "Support for democracy by birth year, no controls", cex=1.2, adj=0.5)
abline(h=seq(-0.4, 0.8, by=0.4), lwd=0.75, lty=3, col=rgb(0.5, 0.5, 0.5, 0.7))
abline(v=seq(1910, 1990, by=10), lwd=0.75, lty=3, col=rgb(0.5, 0.5, 0.5, 0.7))
lines(x=g.b.w.nc.2.coh.eff[[1]]$birthyr, y=g.b.w.nc.2.coh.eff[[1]]$estimate, lwd=2, 
      col=magma(1, begin=0.3, end=0.3))
polygon(x=c(g.b.w.nc.2.coh.eff[[1]]$birthyr, g.b.w.nc.2.coh.eff[[1]]$birthyr[100:1]),
        y=c(g.b.w.nc.2.coh.eff[[1]]$lower, g.b.w.nc.2.coh.eff[[1]]$upper[100:1]), 
        col=magma(1, begin=0.3, end=0.3, alpha=0.5), border=NA, density=NA)
text(x=1998, y=0.88, labels="WVS", pos=2, cex=1)
box()
dev.off()

# year plot

year_labs = c("1995", "1999", "2006", "2011", "2017")

pred.dat.year = data.frame(birthyr=1961, year=c(1995, 1999, 2006, 2011, 2017), agecat4=3)
g.b.w.nc.2.year.pred.eprob = posterior_epred(g.b.w.nc.2.par, newdata=pred.dat.year, ndraws=1000)
g.b.w.nc.2.year.pred.mean = apply(g.b.w.nc.2.year.pred.eprob, 2, mean)
g.b.w.nc.2.year.pred.l95 = apply(g.b.w.nc.2.year.pred.eprob, 2, quantile, probs=0.025)
g.b.w.nc.2.year.pred.u95 = apply(g.b.w.nc.2.year.pred.eprob, 2, quantile, probs=0.975)

pdf("fig4b.pdf", height=4, width=5.5)
par(mfrow=c(1,1), mar=c(3, 3, 2.2, 0.5), tcl=-0.2, las=0, cex=1)
plot(x=1:5, y=g.b.w.nc.2.year.pred.mean, type="n", axes=FALSE, xlab="", ylab="", 
     xlim=c(0.8, 5.2), ylim=c(-0.45, 0.4))
axis(side=1, at=1:5, labels=year_labs, mgp=c(1, 0.1, 0), cex.axis=0.99)
axis(side=2, at=seq(-0.4, 0.4, by=0.2), mgp=c(1, 0.3, 0), cex.axis=0.99)
mtext(side=1, line=1.3, "Year", cex=1.1)
mtext(side=2, line=1.4, "Support democracy & reject autocracy, std scale", cex=1)
mtext(side=3, line=0.5, "Support for democracy by year, no controls", cex=1.2, adj=0.5)
abline(h=seq(-0.4, 0.4, by=0.2), lwd=0.75, lty=3, col=rgb(0.5, 0.5, 0.5, 0.7))
abline(v=1:5, lwd=0.75, lty=3, col=rgb(0.5, 0.5, 0.5, 0.7))
points(x=1:5, y=g.b.w.nc.2.year.pred.mean, pch=16, 
       col=magma(1, begin=0.3, end=0.3), cex=1.2)
segments(x0=1:5, y0=g.b.w.nc.2.year.pred.l95, y1=g.b.w.nc.2.year.pred.u95, col=magma(1, begin=0.3, end=0.3), lwd=1.5)
text(x=5.2, y=0.4, labels="WVS", pos=2, cex=1, adj=1)
box()
dev.off()


## Bayesian GAMs - support democracy & reject autocracy - with controls - in paper

wvs2 = wvs
wvs2$year = as.factor(wvs2$year)
wvs2$agecat4 = as.factor(wvs2$agecat4)
g.b.w.c.2.ser = brm(sup_dem | weights(wt) ~ s(birthyr) + republican + democrat + degree + female + white 
                    + south + income_std + year + agecat4, data=wvs2, prior=g.b.w.2.prior, 
                     control=list(adapt_delta=0.99), chains=6, iter=12, backend="cmdstanr")

g.b.w.c.2.par <- update(
  g.b.w.c.2.ser, chains=3, cores=3, iter=2000, backend="cmdstanr", threads=threading(4),
  control=c(adapt_delta=0.99, max_treedepth=12), stan_model_args=list(stanc_options = list("O1"))
)

summary(g.b.w.c.2.par)
pp_check(g.b.w.c.2.par)
pp_check(g.b.w.c.2.par, type="ecdf_overlay")

pdf("fig_s8b.pdf", width=6, height=6)
plot(g.b.w.c.2.par, newpage=FALSE, 
     variable=c("b_year1999", "b_agecat42", "bs_sbirthyr_1", "sds_sbirthyr_1", "sigma"))
dev.off()

sum.g.b.w.c.2.par <- capture.output(summary(g.b.w.c.2.par))
cat(title="WVS GAM results", sum.g.b.w.c.2.par, file="table3_2.txt", sep="\n", append=TRUE)


# age plot

pred.dat.age = data.frame(birthyr=1961, republican=0, democrat=0, degree=0, female=1, white=1, income_std=0, 
                          south=0, year=2006, agecat4=1:4)
g.b.w.c.2.age.pred.eprob = posterior_epred(g.b.w.c.2.par, newdata=pred.dat.age, ndraws=1000)
g.b.w.c.2.age.pred.mean = apply(g.b.w.c.2.age.pred.eprob, 2, mean)
g.b.w.c.2.age.pred.l95 = apply(g.b.w.c.2.age.pred.eprob, 2, quantile, probs=0.025)
g.b.w.c.2.age.pred.u95 = apply(g.b.w.c.2.age.pred.eprob, 2, quantile, probs=0.975)

pdf("fig3d.pdf", height=4, width=5.5)
par(mfrow=c(1,1), mar=c(3, 3, 2.2, 0.5), tcl=-0.2, las=0, cex=1)
plot(x=1:4, y=g.b.w.c.2.age.pred.mean, type="n", axes=FALSE, xlab="", ylab="", 
     xlim=c(0.8, 4.2), ylim=c(-0.4, 0.2))
axis(side=1, at=1:4, labels=c("18-29", "30-44", "45-59", "60+"), mgp=c(1, 0.1, 0), cex.axis=0.99)
axis(side=2, at=seq(-0.4, 0.2, by=0.2), mgp=c(1, 0.3, 0), cex.axis=0.99)
mtext(side=1, line=1.3, "Age group", cex=1.1)
mtext(side=2, line=1.4, "Support democracy & reject autocracy, std scale", cex=1)
mtext(side=3, line=0.5, "Support for democracy by age group, controls", cex=1.2, adj=0.5)
abline(h=seq(-0.4, 0.2, by=0.2), lwd=0.75, lty=3, col=rgb(0.5, 0.5, 0.5, 0.7))
abline(v=1:4, lwd=0.75, lty=3, col=rgb(0.5, 0.5, 0.5, 0.7))
points(x=1:4, y=g.b.w.c.2.age.pred.mean, pch=16, 
       col=magma(1, begin=0.3, end=0.3), cex=1.2)
segments(x0=1:4, y0=g.b.w.c.2.age.pred.l95, y1=g.b.w.c.2.age.pred.u95, col=magma(1, begin=0.3, end=0.3), lwd=1.5)
text(x=4.2, y=0.2, labels="WVS", pos=2, cex=1, adj=1)
box()
dev.off()

# cohort plot

g.b.w.c.2.coh.eff = conditional_effects(g.b.w.c.2.par, "birthyr")

pdf("fig2d.pdf", height=4, width=5.5)
par(mfrow=c(1,1), mar=c(3, 3, 2.2, 0.5), tcl=-0.2, las=0, cex=1)
plot(x=g.b.w.c.2.coh.eff[[1]]$birthyr, y=g.b.w.c.2.coh.eff[[1]]$estimate, type="n", axes=FALSE, 
     xlab="", ylab="", ylim=c(-0.4, 0.9), xaxs="i")
axis(side=1, at=seq(1900, 2000, by=10), mgp=c(1, 0.1, 0), cex.axis=0.99)
axis(side=2, at=seq(-0.4, 0.8, by=0.4), mgp=c(1, 0.3, 0), cex.axis=0.99)
mtext(side=1, line=1.3, "Year of birth", cex=1.1)
mtext(side=2, line=1.4, "Support democracy & reject autocracy, std scale", cex=1)
mtext(side=3, line=0.5, "Support for democracy by birth year, controls", cex=1.2, adj=0.5)
abline(h=seq(-0.4, 0.8, by=0.4), lwd=0.75, lty=3, col=rgb(0.5, 0.5, 0.5, 0.7))
abline(v=seq(1910, 1990, by=10), lwd=0.75, lty=3, col=rgb(0.5, 0.5, 0.5, 0.7))
lines(x=g.b.w.c.2.coh.eff[[1]]$birthyr, y=g.b.w.c.2.coh.eff[[1]]$estimate, lwd=2, 
      col=magma(1, begin=0.3, end=0.3))
polygon(x=c(g.b.w.c.2.coh.eff[[1]]$birthyr, g.b.w.c.2.coh.eff[[1]]$birthyr[100:1]),
        y=c(g.b.w.c.2.coh.eff[[1]]$lower, g.b.w.c.2.coh.eff[[1]]$upper[100:1]), 
        col=magma(1, begin=0.3, end=0.3, alpha=0.5), border=NA, density=NA)
text(x=1998, y=0.88, labels="WVS", pos=2, cex=1)
box()
dev.off()

# year plot

pred.dat.year = data.frame(birthyr=1961, republican=0, democrat=0, degree=0, female=1, white=1, income_std=0, 
                           south=0, year=c(1995, 1999, 2006, 2011, 2017), agecat4=3)
g.b.w.c.2.year.pred.eprob = posterior_epred(g.b.w.c.2.par, newdata=pred.dat.year, ndraws=1000)
g.b.w.c.2.year.pred.mean = apply(g.b.w.c.2.year.pred.eprob, 2, mean)
g.b.w.c.2.year.pred.l95 = apply(g.b.w.c.2.year.pred.eprob, 2, quantile, probs=0.025)
g.b.w.c.2.year.pred.u95 = apply(g.b.w.c.2.year.pred.eprob, 2, quantile, probs=0.975)

year_labs = c("1995", "1999", "2006", "2011", "2017")

pdf("fig4d.pdf", height=4, width=5.5)
par(mfrow=c(1,1), mar=c(3, 3, 2.2, 0.5), tcl=-0.2, las=0, cex=1)
plot(x=1:5, y=g.b.w.c.2.year.pred.mean, type="n", axes=FALSE, xlab="", ylab="", 
     xlim=c(0.8, 5.2), ylim=c(-0.45, 0.4))
axis(side=1, at=1:5, labels=year_labs, mgp=c(1, 0.1, 0), cex.axis=0.99)
axis(side=2, at=seq(-0.4, 0.4, by=0.2), mgp=c(1, 0.3, 0), cex.axis=0.99)
mtext(side=1, line=1.3, "Year", cex=1.1)
mtext(side=2, line=1.4, "Support democracy & reject autocracy, std scale", cex=1)
mtext(side=3, line=0.5, "Support for democracy by year, controls", cex=1.2, adj=0.5)
abline(h=seq(-0.4, 0.4, by=0.2), lwd=0.75, lty=3, col=rgb(0.5, 0.5, 0.5, 0.7))
abline(v=1:5, lwd=0.75, lty=3, col=rgb(0.5, 0.5, 0.5, 0.7))
points(x=1:5, y=g.b.w.c.2.year.pred.mean, pch=16, 
       col=magma(1, begin=0.3, end=0.3), cex=1.2)
segments(x0=1:5, y0=g.b.w.c.2.year.pred.l95, y1=g.b.w.c.2.year.pred.u95, col=magma(1, begin=0.3, end=0.3), lwd=1.5)
text(x=5.2, y=0.4, labels="WVS", pos=2, cex=1, adj=1)
box()
dev.off()







